Dataset

  • Training HC (N=541)
  • KU HC (N=209)
  • CNUH patient (N=233)
    • TRS : 46
    • FE-SSD : 72
  • Input data : the 77 gray matter FreeSurfer features.
    1. ICV : 1 variable
    2. Subcortical volume : 8 variables
    3. Mean cortical thickness : 34 variables
    4. Mean cortical surface area : 34 variables
      Total 77 variables


Strategy of brain age prediction

  1. Sample exclusion criteria

    • Exclude samples below 18 years and above 75 years old.
    • Exclude samples included in scanning sites with less than 20 control samples.
  2. Data partitioning

    • Train HC Control : Healthy controls from CNUH+Beijing.
    • Test HC Control : Healthy control from KU.
    • Test Patient : Patient from CNUH.
  3. Model development

    • Perform Combat to remove site effect.
    • After removing site effect, build brain age prediction model based on the 77 gray matter FreeSurfer features using Ridge, SVR, RVR, 10-fold Nested cross-validation (10 time).
  4. Test prediction model in Test control data, CNUH patient using variance explained on cross-validation (VEcv), Mean absoulte error (MAE), Peason correlation coefficient (\(r\)).

  5. Calculate Brain-PAD (predicted Age - chronological Age) and corrected Brain-PAD (Goodness-of-fit test)


Data Process


library(openxlsx)
library(dplyr)
## 
## 다음의 패키지를 부착합니다: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
whole_dat <- read.xlsx("2023-1114-BrainAge-sMRI-Final.xlsx",sheet=1,startRow = 1)
Sample exclusion with Age & Site
# Age (filtering subject with age <18 or >75)
whole_dat$Age <- as.numeric(whole_dat$Age)
whole_dat$Age <- floor(whole_dat$Age)
sum(whole_dat$Age < 18 | whole_dat$Age > 75)
## [1] 93
table(whole_dat$Site[which(whole_dat$Age < 18 | whole_dat$Age > 75)])
## 
## IPCAS   SWU 
##    92     1
whole_dat <- whole_dat[-which(whole_dat$Age < 18 | whole_dat$Age > 75),] #1075

# Site
table(whole_dat$Site)
## 
##    BNU   CNUH    HNU  IACAS  IPCAS   JHNU     KU    SWU XHCUMS 
##     56    523      3      2    100     30    214     70     23
whole_dat <- whole_dat[-which(whole_dat$Site=="HNU" | whole_dat$Site=="IACAS"),]

table(whole_dat$Site_2)
## 
##  BEI CNUH   KU 
##  279  523  214
table(whole_dat$Status)
## 
## Control Patient 
##     783     233
# Education
whole_dat$Edu[whole_dat$Edu==1 & !is.na(whole_dat$Edu)] <- 7.5
whole_dat$Edu[whole_dat$Edu==2 & !is.na(whole_dat$Edu)] <- 14
whole_dat$Edu[whole_dat$Edu==3 & !is.na(whole_dat$Edu)] <- 18
whole_dat$Edu <- as.numeric(whole_dat$Edu)
Reduce and seperation Reduce of dataset(KU reslut reference)
table(whole_dat$Status)
## 
## Control Patient 
##     783     233
table(whole_dat$Site_2)
## 
##  BEI CNUH   KU 
##  279  523  214
'%!in%' <- Negate('%in%')
whole_dat <- whole_dat %>% filter(whole_dat$Freesurfer_Number %!in% c("H148","H252","J373"))
whole_dat <- whole_dat %>% filter(whole_dat$Freesurfer_Number %!in% c("K33","K49","K65","K75","K216"))

HC_data <- whole_dat[whole_dat$Status=="Control",]
Patient_data <- setdiff(whole_dat,HC_data)

HC_data <- HC_data[,c(-2)]
Patient_data <- Patient_data[,c(-2)]


Combat


library(sva)
library(factoextra)

whole_dat$Site <- ifelse(whole_dat$Site=="CNUH","JBUH",whole_dat$Site)
whole_dat$Site <- ifelse(whole_dat$Site=="KU","KUAH",whole_dat$Site)

batch <- as.numeric(as.factor(whole_dat$Site))
mod <- model.matrix(~1+Age+Sex+Status,data=whole_dat)
Combat_dat <- ComBat(dat=t(whole_dat[,9:85]), batch=batch, mod=mod)
Combat_dat <- t(Combat_dat)
Combat_dat <- data.frame(SubjID=whole_dat$Freesurfer_Number,Age=whole_dat$Age,
                         Sex=whole_dat$Sex,Site=whole_dat$Site,Combat_dat)

Train_HC_data <- Combat_dat[grep("H|J",whole_dat$Freesurfer_Number),]
Test_HC_data <- Combat_dat[grep("K",whole_dat$Freesurfer_Number),]
Test_Patient_dat <- Combat_dat[grep("C|T",whole_dat$Freesurfer_Number),]

Patient_clinical <- read.xlsx("Patient_clinical_292_0630.xlsx",startRow=1)
Patient_clinical <- Patient_clinical[match(Test_Patient_dat$SubjID,Patient_clinical$freesurfer_ID),]
ID_FESSD <- Patient_clinical$freesurfer_ID[!is.na(Patient_clinical$`FE-SSD`)]
ID_schizo <- Patient_clinical$freesurfer_ID[Patient_clinical$Diagnosis=="schizophrenia"]

Test_TRS_dat <- Test_Patient_dat[grep("T",Test_Patient_dat$SubjID),]
Test_FESSD_dat <- Test_Patient_dat[Test_Patient_dat$SubjID %in% ID_FESSD,]
Test_schizophrenia_dat <- Test_Patient_dat[Test_Patient_dat$SubjID %in% ID_schizo,]

par(mfrow=c(3,2))
plot(density(Train_HC_data$Age,from=15,to=80),main="Age of Train control data")
plot(density(Test_HC_data$Age,from=15,to=80),main="Age of Test control data")
plot(density(Test_Patient_dat$Age,from=15,to=80),main="Age of Test patient data")
plot(density(Test_schizophrenia_dat$Age,from=15,to=80),main="Age of Schizophrenia data")
plot(density(Test_FESSD_dat$Age,from=15,to=80),main="Age of Test FE-SSD data")
plot(density(Test_TRS_dat$Age,from=15,to=80),main="Age of Test TRS data")

PC <- prcomp(scale(whole_dat[,9:85]))
A <- fviz_pca_ind(PC, label="none", habillage=whole_dat$Site,
     addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2",title="Before adjustment") +
  labs(x="Dimension1 (29.0%)", y="Dimension2 (20.7%)")

PC <- prcomp(scale(Combat_dat[,5:81]))
B <- fviz_pca_ind(PC, label="none", habillage=Combat_dat$Site,
     addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2",title="After adjustment") + 
  labs(x="Dimension1 (30.0%)", y="Dimension2 (20.8%)")

library(MASS)
f <- as.formula(paste("Site~",paste(colnames(whole_dat)[9:85],collapse ="+"),sep=""))
lda <- lda(f, data = whole_dat)
lda.value <- predict(lda)
lda.data <- data.frame(Site = whole_dat$Site, lda = lda.value$x)
C <- ggplot(lda.data) + geom_point(aes(lda.LD1, lda.LD2, colour = Site), size = 2.5)

lda <- lda(f, data = Combat_dat)
lda.value <- predict(lda)
lda.data <- data.frame(Site = Combat_dat$Site, lda = lda.value$x)
D <- ggplot(lda.data) + geom_point(aes(lda.LD1, lda.LD2, colour = Site), size = 2.5)
A;B;C;D

library(ggpubr)
png("s1.png", width = 1200, height = 800)
ggarrange(A,B,C,D, 
          labels = c("A","B","C","D"), 
          ncol = 2, nrow = 2,
          font.label = list(size=20))
dev.off()
## png 
##   2
# PCA outlier dection in Train HC data
PC <- prcomp(scale(Train_HC_data[,5:81]))
fviz_pca_ind(PC,col.ind = "#00AFBB",repel=TRUE)
## Warning: ggrepel: 465 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

fviz_pca_ind(PC, addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2", title="") + labs(x="Dimension1 (30.3%)", y="Dimension2 (18.1%)")

out <- c(631,95,171,75,233,646,159,766,223,247,267,5,154,186,192,268,36,238,256,3,249,41,211,197,657) #1067

Train_HC_data$SubjID[rownames(Train_HC_data) %in% out]
##  [1] "H03"  "H05"  "H36"  "H41"  "H75"  "H95"  "H155" "H160" "H172" "H187"
## [11] "H193" "H198" "H212" "H224" "H234" "H239" "H248" "H268" "H250" "H260"
## [21] "H263" "J161" "J248" "J259" "J375"
length(out) # 25
## [1] 25
sum(!rownames(Train_HC_data) %in% out) # 566 --> 541
## [1] 541
Train_HC_data <- Train_HC_data[which(!rownames(Train_HC_data) %in% out),]


Demographics and Clinical characteristics


#Train_HC vs Test_HC vs Schizophrenia
summary_dat_3 <- data.frame(rbind(Train_HC_data,Test_HC_data,Test_schizophrenia_dat))
Edu_dat <- whole_dat[,c(1,6)]
colnames(Edu_dat)[1] <- c("SubjID")
summary_dat_3 <- inner_join(summary_dat_3,Edu_dat,by="SubjID")
summary_dat_3$set <- rep(c(1:3),c(nrow(Train_HC_data),nrow(Test_HC_data),nrow(Test_schizophrenia_dat)))
summary_dat_3$set <- as.factor(summary_dat_3$set)


three.fun_var <- function(var,group,DATA) {
  a11 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==1])),sep="")
  a21 <- paste(round(mean(DATA[var][DATA[group]==1],na.rm = T),2),
               round(sd(DATA[var][DATA[group]==1],na.rm = T),2),sep = " ± ")
  a_1 <- paste(a21,"(",a11,")",sep="")
  
  a12 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==2])),sep="")
  a22 <- paste(round(mean(DATA[var][DATA[group]==2],na.rm = T),2),
               round(sd(DATA[var][DATA[group]==2],na.rm = T),2),sep = " ± ")
  a_2 <-  paste(a22,"(",a12,")",sep="")
  
  a13 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==3])),sep="")
  a23 <- paste(round(mean(DATA[var][DATA[group]==3],na.rm = T),2),
               round(sd(DATA[var][DATA[group]==3],na.rm = T),2),sep = " ± ")
  a_3 <-  paste(a23,"(",a13,")",sep="")
  
  f <- as.formula(paste(var,"~",group))
  Aov <- aov(lm(f,data=DATA))
  b <- round(summary(Aov)[[1]][["Pr(>F)"]][1],4)
  
  y <- DATA[var][,1]
  group <-DATA[group][,1]
  posthoc <- pairwise.t.test(y, group, "bonferroni")
  posthoc <- data.frame(posthoc[3])
  c <- round(posthoc[1,1],4)
  d <- round(posthoc[2,1],4)
  e <- round(posthoc[2,2],4)
  return(matrix(c(a_1,a_2,a_3,b,c,d,e),1,7))
}


TABLE_clinical_3 <- three.fun_var("Age","set",summary_dat_3) 
TABLE_clinical_3 <- rbind(TABLE_clinical_3,three.fun_var("Edu","set",summary_dat_3))
t <- table(summary_dat_3$Sex,summary_dat_3$set)
Sex <- c(paste(t[,1],collapse = "/"),paste(t[,2],collapse = "/"),paste(t[,3],collapse = "/"),round(chisq.test(t)$p.value,4))
Sex <- c(Sex,NA,NA,NA)
TABLE_clinical_3 <- rbind(Sex,TABLE_clinical_3)
rownames(TABLE_clinical_3) <- c("Sex(M/F)","Age","Education, years")
colnames(TABLE_clinical_3) <- c("Health control (CNUH, Beijing)","Health control (KU)","Schizophrenia (CNUH)","p-value","Group 1 vs Group 2", "Group 1 vs Group 3","Group 2 vs Group 3")
TABLE_clinical_3
##                  Health control (CNUH, Beijing) Health control (KU)    
## Sex(M/F)         "281/260"                      "77/132"               
## Age              "31.9 ± 13.7(n=541)"          "37.66 ± 14.66(n=209)"
## Education, years "13.94 ± 2.41(n=267)"         "13.83 ± 2.35(n=209)" 
##                  Schizophrenia (CNUH)   p-value  Group 1 vs Group 2
## Sex(M/F)         "108/88"               "2e-04"  NA                
## Age              "39.34 ± 12.9(n=196)" "0"      "0"               
## Education, years "13.19 ± 3.17(n=195)" "0.0074" "1"               
##                  Group 1 vs Group 3 Group 2 vs Group 3
## Sex(M/F)         NA                 NA                
## Age              "0"                "0.6544"          
## Education, years "0.0084"           "0.0473"
#TestHC vs FE-SSD vs TRS
summary_dat_p3 <- data.frame(rbind(Test_HC_data,Test_FESSD_dat,Test_TRS_dat))
summary_dat_p3 <- inner_join(summary_dat_p3,Edu_dat,by="SubjID")
summary_dat_p3$set <- rep(c(1:3),c(nrow(Test_HC_data),nrow(Test_FESSD_dat),nrow(Test_TRS_dat)))
summary_dat_p3$set <- as.factor(summary_dat_p3$set)


TABLE_clinical_p3 <- three.fun_var("Age","set",summary_dat_p3) 
TABLE_clinical_p3 <- rbind(TABLE_clinical_p3,three.fun_var("Edu","set",summary_dat_p3))
t <- table(summary_dat_p3$Sex,summary_dat_p3$set)
Sex <- c(paste(t[,1],collapse = "/"),paste(t[,2],collapse = "/"),paste(t[,3],collapse = "/"),round(chisq.test(t)$p.value,3))
Sex <- c(Sex,NA,NA,NA)
TABLE_clinical_p3 <- rbind(Sex,TABLE_clinical_p3)
rownames(TABLE_clinical_p3) <- c("Sex(M/F)","Age","Education, years")
colnames(TABLE_clinical_p3) <- c("Test HC","FE-SSD","TRS","p-value","Group 1 vs Group 2", "Group 1 vs Group 3","Group 2 vs Group 3")
TABLE_clinical_p3
##                  Test HC                 FE-SSD                
## Sex(M/F)         "77/132"                "28/44"               
## Age              "37.66 ± 14.66(n=209)" "29.61 ± 12.26(n=72)"
## Education, years "13.83 ± 2.35(n=209)"  "13.26 ± 2.79(n=72)" 
##                  TRS                   p-value  Group 1 vs Group 2
## Sex(M/F)         "30/16"               "0.002"  NA                
## Age              "42.61 ± 9.9(n=46)"  "0"      "1e-04"           
## Education, years "13.73 ± 2.36(n=46)" "0.2405" "0.277"           
##                  Group 1 vs Group 3 Group 2 vs Group 3
## Sex(M/F)         NA                 NA                
## Age              "0.0776"           "0"               
## Education, years "1"                "0.9508"


Comparison of the 77 gray matter Freesurfer features


three.fun_var <- function(var,group,DATA) {
  a11 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==1])),sep="")
  a21 <- paste(round(mean(DATA[var][DATA[group]==1],na.rm = T),2),
               round(sd(DATA[var][DATA[group]==1],na.rm = T),2),sep = " ± ")
  a_1 <- paste(a21,"(",a11,")",sep="")
  
  a12 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==2])),sep="")
  a22 <- paste(round(mean(DATA[var][DATA[group]==2],na.rm = T),2),
               round(sd(DATA[var][DATA[group]==2],na.rm = T),2),sep = " ± ")
  a_2 <-  paste(a22,"(",a12,")",sep="")
  
  a13 <- paste("n=",sum(!is.na(DATA[var][DATA[group]==3])),sep="")
  a23 <- paste(round(mean(DATA[var][DATA[group]==3],na.rm = T),2),
               round(sd(DATA[var][DATA[group]==3],na.rm = T),2),sep = " ± ")
  a_3 <-  paste(a23,"(",a13,")",sep="")
  
  f <- as.formula(paste(var,"~","set"))
  Aov <- aov(lm(f,data=DATA))
  b <- round(summary(Aov)[[1]][["Pr(>F)"]][1],4)

  y <- DATA[var][,1]
  group <-DATA[group][,1]
  posthoc <- pairwise.t.test(y, group, "bonferroni")
  posthoc <- data.frame(posthoc[3])
  c <- round(posthoc[1,1],4)
  d <- round(posthoc[2,1],4)
  e <- round(posthoc[2,2],4)
  
  f1 <- as.formula(paste(var,"~","Age + Sex + Edu + set"))
  Aov1 <- aov(lm(f1,data=DATA))
  f <- round(summary(Aov1)[[1]][["Pr(>F)"]][4],4)
  return(matrix(c(a_1,a_2,a_3,b,c,d,e,f),1,8))
}

conti_var <- setdiff(colnames(summary_dat_3),c("SubjID","Age","Site","Sex","set","Edu"))

TABLE_sMRI_3group <- NULL
for (i in 1:length(conti_var)){
  TABLE_sMRI_3group  <- rbind(TABLE_sMRI_3group,three.fun_var(conti_var[i],"set",summary_dat_3))
}
rownames(TABLE_sMRI_3group) <- conti_var
colnames(TABLE_sMRI_3group) <- c("Health control (CNUH, Beijing)","Health control (KU)","Schizophrenia (CNUH)","p-value","Group 1 vs Group 2", "Group 1 vs Group 3","Group 2 vs Group 3","age.sex.edu.adj-p")
TABLE_sMRI_3group <- data.frame(TABLE_sMRI_3group)
TABLE_sMRI_3group$p.value <- as.numeric(paste(TABLE_sMRI_3group$p.value))
TABLE_sMRI_3group$bonf.adj.p.value <- p.adjust(TABLE_sMRI_3group$p.value,method="bonferroni")
TABLE_sMRI_3group$bonf.age.sex.edu.adj.p.value <- p.adjust(TABLE_sMRI_3group$age.sex.edu.adj.p,method="bonferroni")
TABLE_sMRI_3group
##                                       Health.control..CNUH..Beijing.
## Avg_bankssts_thickavg                            2.62 ± 0.13(n=541)
## Avg_caudalanteriorcingulate_thickavg             2.71 ± 0.18(n=541)
## Avg_caudalmiddlefrontal_thickavg                 2.66 ± 0.11(n=541)
## Avg_cuneus_thickavg                               1.9 ± 0.11(n=541)
## Avg_entorhinal_thickavg                          3.47 ± 0.25(n=541)
## Avg_fusiform_thickavg                             2.8 ± 0.11(n=541)
## Avg_inferiorparietal_thickavg                    2.57 ± 0.11(n=541)
## Avg_inferiortemporal_thickavg                    2.87 ± 0.11(n=541)
## Avg_isthmuscingulate_thickavg                    2.39 ± 0.16(n=541)
## Avg_lateraloccipital_thickavg                      2.3 ± 0.1(n=541)
## Avg_lateralorbitofrontal_thickavg                2.74 ± 0.12(n=541)
## Avg_lingual_thickavg                             2.01 ± 0.11(n=541)
## Avg_medialorbitofrontal_thickavg                 2.58 ± 0.13(n=541)
## Avg_middletemporal_thickavg                      2.97 ± 0.12(n=541)
## Avg_parahippocampal_thickavg                     2.69 ± 0.21(n=541)
## Avg_paracentral_thickavg                         2.51 ± 0.13(n=541)
## Avg_parsopercularis_thickavg                     2.68 ± 0.13(n=541)
## Avg_parsorbitalis_thickavg                        2.8 ± 0.15(n=541)
## Avg_parstriangularis_thickavg                    2.58 ± 0.12(n=541)
## Avg_pericalcarine_thickavg                       1.61 ± 0.11(n=541)
## Avg_postcentral_thickavg                         2.12 ± 0.11(n=541)
## Avg_posteriorcingulate_thickavg                  2.57 ± 0.12(n=541)
## Avg_precentral_thickavg                          2.62 ± 0.12(n=541)
## Avg_precuneus_thickavg                           2.45 ± 0.11(n=541)
## Avg_rostralanteriorcingulate_thickavg            2.95 ± 0.16(n=541)
## Avg_rostralmiddlefrontal_thickavg                 2.5 ± 0.11(n=541)
## Avg_superiorfrontal_thickavg                     2.88 ± 0.12(n=541)
## Avg_superiorparietal_thickavg                     2.28 ± 0.1(n=541)
## Avg_superiortemporal_thickavg                    2.88 ± 0.14(n=541)
## Avg_supramarginal_thickavg                       2.62 ± 0.12(n=541)
## Avg_frontalpole_thickavg                         2.97 ± 0.24(n=541)
## Avg_temporalpole_thickavg                        3.75 ± 0.22(n=541)
## Avg_transversetemporal_thickavg                  2.34 ± 0.17(n=541)
## Avg_insula_thickavg                              3.04 ± 0.14(n=541)
## Avg_bankssts_surfavg                         967.57 ± 121.79(n=541)
## Avg_caudalanteriorcingulate_surfavg          645.48 ± 108.22(n=541)
## Avg_caudalmiddlefrontal_surfavg              2165.97 ± 304.8(n=541)
## Avg_cuneus_surfavg                          1625.94 ± 191.04(n=541)
## Avg_entorhinal_surfavg                        419.46 ± 68.21(n=541)
## Avg_fusiform_surfavg                        3129.33 ± 324.28(n=541)
## Avg_inferiorparietal_surfavg                 5216.56 ± 635.7(n=541)
## Avg_inferiortemporal_surfavg                3494.81 ± 429.24(n=541)
## Avg_isthmuscingulate_surfavg                1008.34 ± 133.08(n=541)
## Avg_lateraloccipital_surfavg                5117.72 ± 579.53(n=541)
## Avg_lateralorbitofrontal_surfavg            2588.52 ± 250.91(n=541)
## Avg_lingual_surfavg                         3231.03 ± 386.93(n=541)
## Avg_medialorbitofrontal_surfavg             1914.23 ± 188.57(n=541)
## Avg_middletemporal_surfavg                   3476.88 ± 401.7(n=541)
## Avg_parahippocampal_surfavg                    662.46 ± 69.9(n=541)
## Avg_paracentral_surfavg                      1432.7 ± 147.16(n=541)
## Avg_parsopercularis_surfavg                 1504.98 ± 196.49(n=541)
## Avg_parsorbitalis_surfavg                      728.13 ± 80.8(n=541)
## Avg_parstriangularis_surfavg                1450.58 ± 191.35(n=541)
## Avg_pericalcarine_surfavg                   1565.86 ± 244.81(n=541)
## Avg_postcentral_surfavg                     4197.01 ± 426.36(n=541)
## Avg_posteriorcingulate_surfavg              1202.75 ± 150.23(n=541)
## Avg_precentral_surfavg                      4980.56 ± 461.28(n=541)
## Avg_precuneus_surfavg                       4073.36 ± 461.12(n=541)
## Avg_rostralanteriorcingulate_surfavg         707.86 ± 109.49(n=541)
## Avg_rostralmiddlefrontal_surfavg            5662.36 ± 719.79(n=541)
## Avg_superiorfrontal_surfavg                 6853.75 ± 716.93(n=541)
## Avg_superiorparietal_surfavg                5494.47 ± 593.85(n=541)
## Avg_superiortemporal_surfavg                3877.44 ± 387.67(n=541)
## Avg_supramarginal_surfavg                   3952.15 ± 525.93(n=541)
## Avg_frontalpole_surfavg                       277.22 ± 26.03(n=541)
## Avg_temporalpole_surfavg                      473.96 ± 49.44(n=541)
## Avg_transversetemporal_surfavg                412.78 ± 52.03(n=541)
## Avg_insula_surfavg                           2458.61 ± 232.1(n=541)
## ICV                                    1510324.43 ± 157050.4(n=541)
## Avg_LatVent                                7516.06 ± 4127.39(n=541)
## Avg_thal                                    8073.94 ± 842.81(n=541)
## Avg_caud                                    3554.28 ± 456.93(n=541)
## Avg_put                                      5148.74 ± 589.6(n=541)
## Avg_pal                                     2020.69 ± 230.34(n=541)
## Avg_hippo                                   4268.04 ± 370.26(n=541)
## Avg_amyg                                    1758.63 ± 188.75(n=541)
## Avg_accumb                                    484.49 ± 81.47(n=541)
##                                                  Health.control..KU.
## Avg_bankssts_thickavg                             2.6 ± 0.13(n=209)
## Avg_caudalanteriorcingulate_thickavg             2.68 ± 0.18(n=209)
## Avg_caudalmiddlefrontal_thickavg                 2.64 ± 0.13(n=209)
## Avg_cuneus_thickavg                              1.88 ± 0.11(n=209)
## Avg_entorhinal_thickavg                          3.46 ± 0.26(n=209)
## Avg_fusiform_thickavg                            2.78 ± 0.13(n=209)
## Avg_inferiorparietal_thickavg                    2.55 ± 0.11(n=209)
## Avg_inferiortemporal_thickavg                    2.85 ± 0.12(n=209)
## Avg_isthmuscingulate_thickavg                    2.36 ± 0.17(n=209)
## Avg_lateraloccipital_thickavg                    2.29 ± 0.11(n=209)
## Avg_lateralorbitofrontal_thickavg                2.71 ± 0.15(n=209)
## Avg_lingual_thickavg                             1.99 ± 0.11(n=209)
## Avg_medialorbitofrontal_thickavg                 2.55 ± 0.15(n=209)
## Avg_middletemporal_thickavg                      2.95 ± 0.13(n=209)
## Avg_parahippocampal_thickavg                     2.68 ± 0.22(n=209)
## Avg_paracentral_thickavg                         2.49 ± 0.14(n=209)
## Avg_parsopercularis_thickavg                     2.65 ± 0.13(n=209)
## Avg_parsorbitalis_thickavg                       2.77 ± 0.18(n=209)
## Avg_parstriangularis_thickavg                    2.55 ± 0.14(n=209)
## Avg_pericalcarine_thickavg                        1.6 ± 0.11(n=209)
## Avg_postcentral_thickavg                          2.1 ± 0.11(n=209)
## Avg_posteriorcingulate_thickavg                  2.55 ± 0.13(n=209)
## Avg_precentral_thickavg                           2.6 ± 0.13(n=209)
## Avg_precuneus_thickavg                           2.43 ± 0.12(n=209)
## Avg_rostralanteriorcingulate_thickavg            2.92 ± 0.17(n=209)
## Avg_rostralmiddlefrontal_thickavg                2.48 ± 0.12(n=209)
## Avg_superiorfrontal_thickavg                     2.86 ± 0.14(n=209)
## Avg_superiorparietal_thickavg                    2.26 ± 0.11(n=209)
## Avg_superiortemporal_thickavg                    2.85 ± 0.14(n=209)
## Avg_supramarginal_thickavg                        2.6 ± 0.13(n=209)
## Avg_frontalpole_thickavg                         2.95 ± 0.24(n=209)
## Avg_temporalpole_thickavg                        3.73 ± 0.22(n=209)
## Avg_transversetemporal_thickavg                   2.3 ± 0.19(n=209)
## Avg_insula_thickavg                                 3 ± 0.16(n=209)
## Avg_bankssts_surfavg                         944.64 ± 135.77(n=209)
## Avg_caudalanteriorcingulate_surfavg          631.71 ± 121.04(n=209)
## Avg_caudalmiddlefrontal_surfavg             2108.21 ± 324.36(n=209)
## Avg_cuneus_surfavg                          1595.19 ± 209.57(n=209)
## Avg_entorhinal_surfavg                         417.4 ± 68.72(n=209)
## Avg_fusiform_surfavg                           3070.39 ± 380(n=209)
## Avg_inferiorparietal_surfavg                 5090.6 ± 691.75(n=209)
## Avg_inferiortemporal_surfavg                3422.45 ± 477.41(n=209)
## Avg_isthmuscingulate_surfavg                 985.27 ± 146.75(n=209)
## Avg_lateraloccipital_surfavg                 5010.8 ± 641.02(n=209)
## Avg_lateralorbitofrontal_surfavg            2545.06 ± 280.94(n=209)
## Avg_lingual_surfavg                         3176.06 ± 405.42(n=209)
## Avg_medialorbitofrontal_surfavg             1889.05 ± 214.34(n=209)
## Avg_middletemporal_surfavg                  3389.74 ± 457.63(n=209)
## Avg_parahippocampal_surfavg                    652.13 ± 78.2(n=209)
## Avg_paracentral_surfavg                     1418.23 ± 164.58(n=209)
## Avg_parsopercularis_surfavg                 1474.17 ± 207.84(n=209)
## Avg_parsorbitalis_surfavg                      713.9 ± 89.72(n=209)
## Avg_parstriangularis_surfavg                1414.59 ± 212.28(n=209)
## Avg_pericalcarine_surfavg                   1544.01 ± 254.77(n=209)
## Avg_postcentral_surfavg                     4132.66 ± 475.47(n=209)
## Avg_posteriorcingulate_surfavg              1175.53 ± 170.86(n=209)
## Avg_precentral_surfavg                      4904.58 ± 519.49(n=209)
## Avg_precuneus_surfavg                       3983.16 ± 522.64(n=209)
## Avg_rostralanteriorcingulate_surfavg         691.98 ± 130.93(n=209)
## Avg_rostralmiddlefrontal_surfavg             5501.88 ± 839.1(n=209)
## Avg_superiorfrontal_surfavg                 6706.45 ± 847.65(n=209)
## Avg_superiorparietal_surfavg                5380.33 ± 666.04(n=209)
## Avg_superiortemporal_surfavg                3799.53 ± 458.77(n=209)
## Avg_supramarginal_surfavg                   3868.51 ± 590.94(n=209)
## Avg_frontalpole_surfavg                       272.64 ± 29.73(n=209)
## Avg_temporalpole_surfavg                      470.19 ± 50.68(n=209)
## Avg_transversetemporal_surfavg                405.54 ± 56.68(n=209)
## Avg_insula_surfavg                          2427.76 ± 267.23(n=209)
## ICV                                   1479141.93 ± 178127.94(n=209)
## Avg_LatVent                                7859.21 ± 4663.11(n=209)
## Avg_thal                                    7848.64 ± 875.31(n=209)
## Avg_caud                                       3460 ± 458.02(n=209)
## Avg_put                                     4990.17 ± 616.09(n=209)
## Avg_pal                                      1983.01 ± 233.1(n=209)
## Avg_hippo                                   4188.57 ± 399.62(n=209)
## Avg_amyg                                     1720.7 ± 204.65(n=209)
## Avg_accumb                                     469.26 ± 84.9(n=209)
##                                                 Schizophrenia..CNUH. p.value
## Avg_bankssts_thickavg                            2.52 ± 0.12(n=196)  0.0000
## Avg_caudalanteriorcingulate_thickavg             2.65 ± 0.19(n=196)  0.0016
## Avg_caudalmiddlefrontal_thickavg                 2.54 ± 0.12(n=196)  0.0000
## Avg_cuneus_thickavg                              1.86 ± 0.11(n=196)  0.0000
## Avg_entorhinal_thickavg                           3.4 ± 0.25(n=196)  0.0106
## Avg_fusiform_thickavg                            2.71 ± 0.12(n=196)  0.0000
## Avg_inferiorparietal_thickavg                    2.47 ± 0.11(n=196)  0.0000
## Avg_inferiortemporal_thickavg                    2.78 ± 0.12(n=196)  0.0000
## Avg_isthmuscingulate_thickavg                    2.27 ± 0.16(n=196)  0.0000
## Avg_lateraloccipital_thickavg                     2.22 ± 0.1(n=196)  0.0000
## Avg_lateralorbitofrontal_thickavg                 2.6 ± 0.13(n=196)  0.0000
## Avg_lingual_thickavg                              1.98 ± 0.1(n=196)  0.0011
## Avg_medialorbitofrontal_thickavg                 2.47 ± 0.14(n=196)  0.0000
## Avg_middletemporal_thickavg                      2.86 ± 0.13(n=196)  0.0000
## Avg_parahippocampal_thickavg                      2.6 ± 0.23(n=196)  0.0000
## Avg_paracentral_thickavg                         2.44 ± 0.12(n=196)  0.0000
## Avg_parsopercularis_thickavg                     2.55 ± 0.14(n=196)  0.0000
## Avg_parsorbitalis_thickavg                       2.64 ± 0.16(n=196)  0.0000
## Avg_parstriangularis_thickavg                    2.45 ± 0.12(n=196)  0.0000
## Avg_pericalcarine_thickavg                        1.58 ± 0.1(n=196)  0.0027
## Avg_postcentral_thickavg                         2.03 ± 0.09(n=196)  0.0000
## Avg_posteriorcingulate_thickavg                  2.48 ± 0.14(n=196)  0.0000
## Avg_precentral_thickavg                          2.52 ± 0.13(n=196)  0.0000
## Avg_precuneus_thickavg                           2.36 ± 0.11(n=196)  0.0000
## Avg_rostralanteriorcingulate_thickavg            2.84 ± 0.18(n=196)  0.0000
## Avg_rostralmiddlefrontal_thickavg                 2.37 ± 0.1(n=196)  0.0000
## Avg_superiorfrontal_thickavg                     2.76 ± 0.12(n=196)  0.0000
## Avg_superiorparietal_thickavg                      2.2 ± 0.1(n=196)  0.0000
## Avg_superiortemporal_thickavg                    2.74 ± 0.14(n=196)  0.0000
## Avg_supramarginal_thickavg                       2.52 ± 0.12(n=196)  0.0000
## Avg_frontalpole_thickavg                         2.79 ± 0.23(n=196)  0.0000
## Avg_temporalpole_thickavg                        3.62 ± 0.24(n=196)  0.0000
## Avg_transversetemporal_thickavg                  2.28 ± 0.17(n=196)  0.0002
## Avg_insula_thickavg                              2.94 ± 0.14(n=196)  0.0000
## Avg_bankssts_surfavg                          948.13 ± 137.1(n=196)  0.0406
## Avg_caudalanteriorcingulate_surfavg          597.62 ± 107.91(n=196)  0.0000
## Avg_caudalmiddlefrontal_surfavg             2105.85 ± 317.56(n=196)  0.0163
## Avg_cuneus_surfavg                          1637.13 ± 211.18(n=196)  0.0784
## Avg_entorhinal_surfavg                        417.58 ± 66.72(n=196)  0.9069
## Avg_fusiform_surfavg                        3001.54 ± 360.22(n=196)  0.0000
## Avg_inferiorparietal_surfavg                5120.73 ± 659.51(n=196)  0.0314
## Avg_inferiortemporal_surfavg                3437.63 ± 457.53(n=196)  0.0812
## Avg_isthmuscingulate_surfavg                1011.05 ± 139.99(n=196)  0.0860
## Avg_lateraloccipital_surfavg                5050.99 ± 633.71(n=196)  0.0716
## Avg_lateralorbitofrontal_surfavg            2600.96 ± 291.54(n=196)  0.0703
## Avg_lingual_surfavg                         3170.42 ± 396.38(n=196)  0.0821
## Avg_medialorbitofrontal_surfavg             1902.11 ± 215.58(n=196)  0.2893
## Avg_middletemporal_surfavg                  3416.84 ± 442.17(n=196)  0.0238
## Avg_parahippocampal_surfavg                   642.94 ± 75.34(n=196)  0.0040
## Avg_paracentral_surfavg                     1443.73 ± 171.59(n=196)  0.2557
## Avg_parsopercularis_surfavg                 1452.61 ± 223.15(n=196)  0.0054
## Avg_parsorbitalis_surfavg                     733.48 ± 96.98(n=196)  0.0528
## Avg_parstriangularis_surfavg                1425.49 ± 224.21(n=196)  0.0619
## Avg_pericalcarine_surfavg                   1553.84 ± 242.96(n=196)  0.5293
## Avg_postcentral_surfavg                        4163 ± 469.17(n=196)  0.1896
## Avg_posteriorcingulate_surfavg              1167.49 ± 169.67(n=196)  0.0107
## Avg_precentral_surfavg                      4967.66 ± 519.56(n=196)  0.1567
## Avg_precuneus_surfavg                       4018.58 ± 496.78(n=196)  0.0542
## Avg_rostralanteriorcingulate_surfavg          678.14 ± 129.5(n=196)  0.0079
## Avg_rostralmiddlefrontal_surfavg             5744.03 ± 864.2(n=196)  0.0054
## Avg_superiorfrontal_surfavg                 6819.75 ± 835.58(n=196)  0.0648
## Avg_superiorparietal_surfavg                5442.31 ± 631.79(n=196)  0.0710
## Avg_superiortemporal_surfavg                3800.74 ± 431.61(n=196)  0.0177
## Avg_supramarginal_surfavg                   3906.44 ± 564.67(n=196)  0.1533
## Avg_frontalpole_surfavg                       276.63 ± 31.96(n=196)  0.1308
## Avg_temporalpole_surfavg                       479.5 ± 54.55(n=196)  0.1792
## Avg_transversetemporal_surfavg                394.04 ± 54.84(n=196)  0.0001
## Avg_insula_surfavg                          2416.67 ± 262.66(n=196)  0.0750
## ICV                                   1508777.92 ± 155869.81(n=196)  0.0526
## Avg_LatVent                                   9855 ± 5140.15(n=196)  0.0000
## Avg_thal                                    7204.15 ± 879.25(n=196)  0.0000
## Avg_caud                                    3391.64 ± 425.19(n=196)  0.0000
## Avg_put                                     4971.56 ± 637.67(n=196)  0.0001
## Avg_pal                                     2096.43 ± 234.23(n=196)  0.0000
## Avg_hippo                                   3945.11 ± 435.22(n=196)  0.0000
## Avg_amyg                                    1664.91 ± 230.46(n=196)  0.0000
## Avg_accumb                                    455.59 ± 84.96(n=196)  0.0001
##                                       Group.1.vs.Group.2 Group.1.vs.Group.3
## Avg_bankssts_thickavg                             0.1927                  0
## Avg_caudalanteriorcingulate_thickavg              0.4142             0.0012
## Avg_caudalmiddlefrontal_thickavg                  0.0861                  0
## Avg_cuneus_thickavg                               0.1531                  0
## Avg_entorhinal_thickavg                                1             0.0082
## Avg_fusiform_thickavg                             0.0721                  0
## Avg_inferiorparietal_thickavg                     0.0456                  0
## Avg_inferiortemporal_thickavg                     0.3127                  0
## Avg_isthmuscingulate_thickavg                     0.2106                  0
## Avg_lateraloccipital_thickavg                     0.3895                  0
## Avg_lateralorbitofrontal_thickavg                 0.0178                  0
## Avg_lingual_thickavg                              0.0761             0.0018
## Avg_medialorbitofrontal_thickavg                   0.036                  0
## Avg_middletemporal_thickavg                         0.08                  0
## Avg_parahippocampal_thickavg                           1                  0
## Avg_paracentral_thickavg                          0.2343                  0
## Avg_parsopercularis_thickavg                      0.0115                  0
## Avg_parsorbitalis_thickavg                        0.0928                  0
## Avg_parstriangularis_thickavg                     0.0363                  0
## Avg_pericalcarine_thickavg                        0.2684             0.0025
## Avg_postcentral_thickavg                           0.255                  0
## Avg_posteriorcingulate_thickavg                   0.0478                  0
## Avg_precentral_thickavg                           0.1142                  0
## Avg_precuneus_thickavg                            0.0431                  0
## Avg_rostralanteriorcingulate_thickavg              0.149                  0
## Avg_rostralmiddlefrontal_thickavg                 0.1316                  0
## Avg_superiorfrontal_thickavg                      0.0561                  0
## Avg_superiorparietal_thickavg                     0.2885                  0
## Avg_superiortemporal_thickavg                     0.0112                  0
## Avg_supramarginal_thickavg                        0.0621                  0
## Avg_frontalpole_thickavg                          0.5615                  0
## Avg_temporalpole_thickavg                         0.8051                  0
## Avg_transversetemporal_thickavg                   0.0783              2e-04
## Avg_insula_thickavg                               0.0017                  0
## Avg_bankssts_surfavg                              0.0851              0.208
## Avg_caudalanteriorcingulate_surfavg               0.3857                  0
## Avg_caudalmiddlefrontal_surfavg                   0.0696             0.0629
## Avg_cuneus_surfavg                                0.1763                  1
## Avg_entorhinal_surfavg                                 1                  1
## Avg_fusiform_surfavg                              0.1083                  0
## Avg_inferiorparietal_surfavg                      0.0544             0.2365
## Avg_inferiortemporal_surfavg                      0.1402             0.3737
## Avg_isthmuscingulate_surfavg                      0.1196                  1
## Avg_lateraloccipital_surfavg                      0.0908             0.5583
## Avg_lateralorbitofrontal_surfavg                  0.1367                  1
## Avg_lingual_surfavg                               0.2587              0.194
## Avg_medialorbitofrontal_surfavg                   0.3689                  1
## Avg_middletemporal_surfavg                        0.0348             0.2672
## Avg_parahippocampal_surfavg                       0.2468             0.0041
## Avg_paracentral_surfavg                           0.7685                  1
## Avg_parsopercularis_surfavg                        0.195             0.0067
## Avg_parsorbitalis_surfavg                           0.13                  1
## Avg_parstriangularis_surfavg                      0.0898             0.4171
## Avg_pericalcarine_surfavg                         0.8308                  1
## Avg_postcentral_surfavg                           0.2316                  1
## Avg_posteriorcingulate_surfavg                    0.1078              0.024
## Avg_precentral_surfavg                            0.1672                  1
## Avg_precuneus_surfavg                              0.066             0.5214
## Avg_rostralanteriorcingulate_surfavg              0.3032             0.0083
## Avg_rostralmiddlefrontal_surfavg                  0.0347             0.6262
## Avg_superiorfrontal_surfavg                       0.0584                  1
## Avg_superiorparietal_surfavg                      0.0709             0.9356
## Avg_superiortemporal_surfavg                      0.0628              0.079
## Avg_supramarginal_surfavg                         0.1851             0.9545
## Avg_frontalpole_surfavg                           0.1386                  1
## Avg_temporalpole_surfavg                               1             0.5736
## Avg_transversetemporal_surfavg                    0.2938              1e-04
## Avg_insula_surfavg                                0.3751             0.1252
## ICV                                               0.0543                  1
## Avg_LatVent                                            1                  0
## Avg_thal                                          0.0039                  0
## Avg_caud                                          0.0312              1e-04
## Avg_put                                           0.0041             0.0014
## Avg_pal                                           0.1387              3e-04
## Avg_hippo                                         0.0383                  0
## Avg_amyg                                          0.0632                  0
## Avg_accumb                                        0.0732              1e-04
##                                       Group.2.vs.Group.3 age.sex.edu.adj.p
## Avg_bankssts_thickavg                                  0                 0
## Avg_caudalanteriorcingulate_thickavg              0.2314            0.0394
## Avg_caudalmiddlefrontal_thickavg                       0                 0
## Avg_cuneus_thickavg                               0.0326            0.0027
## Avg_entorhinal_thickavg                           0.1047            0.0086
## Avg_fusiform_thickavg                                  0                 0
## Avg_inferiorparietal_thickavg                          0                 0
## Avg_inferiortemporal_thickavg                          0                 0
## Avg_isthmuscingulate_thickavg                          0                 0
## Avg_lateraloccipital_thickavg                          0                 0
## Avg_lateralorbitofrontal_thickavg                      0                 0
## Avg_lingual_thickavg                              0.8757            0.1268
## Avg_medialorbitofrontal_thickavg                       0                 0
## Avg_middletemporal_thickavg                            0                 0
## Avg_parahippocampal_thickavg                      0.0015            0.0012
## Avg_paracentral_thickavg                           1e-04                 0
## Avg_parsopercularis_thickavg                           0                 0
## Avg_parsorbitalis_thickavg                             0                 0
## Avg_parstriangularis_thickavg                          0                 0
## Avg_pericalcarine_thickavg                        0.4702            0.0388
## Avg_postcentral_thickavg                               0                 0
## Avg_posteriorcingulate_thickavg                        0                 0
## Avg_precentral_thickavg                                0                 0
## Avg_precuneus_thickavg                                 0                 0
## Avg_rostralanteriorcingulate_thickavg                  0                 0
## Avg_rostralmiddlefrontal_thickavg                      0                 0
## Avg_superiorfrontal_thickavg                           0                 0
## Avg_superiorparietal_thickavg                          0                 0
## Avg_superiortemporal_thickavg                          0                 0
## Avg_supramarginal_thickavg                             0                 0
## Avg_frontalpole_thickavg                               0                 0
## Avg_temporalpole_thickavg                              0                 0
## Avg_transversetemporal_thickavg                   0.3634            0.1065
## Avg_insula_thickavg                                    0                 0
## Avg_bankssts_surfavg                                   1            0.7081
## Avg_caudalanteriorcingulate_surfavg               0.0063                 0
## Avg_caudalmiddlefrontal_surfavg                        1            0.5272
## Avg_cuneus_surfavg                                0.1043            0.3277
## Avg_entorhinal_surfavg                                 1            0.6574
## Avg_fusiform_surfavg                              0.1348             9e-04
## Avg_inferiorparietal_surfavg                           1            0.7843
## Avg_inferiortemporal_surfavg                           1             0.461
## Avg_isthmuscingulate_surfavg                      0.1799             0.441
## Avg_lateraloccipital_surfavg                           1            0.7346
## Avg_lateralorbitofrontal_surfavg                  0.1056            0.1934
## Avg_lingual_surfavg                                    1            0.5824
## Avg_medialorbitofrontal_surfavg                        1            0.4816
## Avg_middletemporal_surfavg                             1            0.6729
## Avg_parahippocampal_surfavg                       0.6169            0.0469
## Avg_paracentral_surfavg                           0.3043            0.6717
## Avg_parsopercularis_surfavg                       0.8696            0.1417
## Avg_parsorbitalis_surfavg                         0.0684            0.2962
## Avg_parstriangularis_surfavg                           1            0.8697
## Avg_pericalcarine_surfavg                              1            0.9883
## Avg_postcentral_surfavg                                1            0.5574
## Avg_posteriorcingulate_surfavg                         1            0.2626
## Avg_precentral_surfavg                            0.5791            0.7428
## Avg_precuneus_surfavg                                  1            0.8551
## Avg_rostralanteriorcingulate_surfavg               0.725            0.0086
## Avg_rostralmiddlefrontal_surfavg                  0.0055             0.013
## Avg_superiorfrontal_surfavg                        0.422            0.5975
## Avg_superiorparietal_surfavg                       0.941            0.7149
## Avg_superiortemporal_surfavg                           1             0.337
## Avg_supramarginal_surfavg                              1            0.5834
## Avg_frontalpole_surfavg                           0.4631             0.756
## Avg_temporalpole_surfavg                           0.197            0.5612
## Avg_transversetemporal_surfavg                    0.0946            0.0021
## Avg_insula_surfavg                                     1            0.0757
## ICV                                               0.1968            0.7287
## Avg_LatVent                                            0             2e-04
## Avg_thal                                               0                 0
## Avg_caud                                          0.3826            0.0188
## Avg_put                                                1            0.1725
## Avg_pal                                                0                 0
## Avg_hippo                                              0                 0
## Avg_amyg                                          0.0164                 0
## Avg_accumb                                        0.2939            0.0614
##                                       bonf.adj.p.value
## Avg_bankssts_thickavg                           0.0000
## Avg_caudalanteriorcingulate_thickavg            0.1232
## Avg_caudalmiddlefrontal_thickavg                0.0000
## Avg_cuneus_thickavg                             0.0000
## Avg_entorhinal_thickavg                         0.8162
## Avg_fusiform_thickavg                           0.0000
## Avg_inferiorparietal_thickavg                   0.0000
## Avg_inferiortemporal_thickavg                   0.0000
## Avg_isthmuscingulate_thickavg                   0.0000
## Avg_lateraloccipital_thickavg                   0.0000
## Avg_lateralorbitofrontal_thickavg               0.0000
## Avg_lingual_thickavg                            0.0847
## Avg_medialorbitofrontal_thickavg                0.0000
## Avg_middletemporal_thickavg                     0.0000
## Avg_parahippocampal_thickavg                    0.0000
## Avg_paracentral_thickavg                        0.0000
## Avg_parsopercularis_thickavg                    0.0000
## Avg_parsorbitalis_thickavg                      0.0000
## Avg_parstriangularis_thickavg                   0.0000
## Avg_pericalcarine_thickavg                      0.2079
## Avg_postcentral_thickavg                        0.0000
## Avg_posteriorcingulate_thickavg                 0.0000
## Avg_precentral_thickavg                         0.0000
## Avg_precuneus_thickavg                          0.0000
## Avg_rostralanteriorcingulate_thickavg           0.0000
## Avg_rostralmiddlefrontal_thickavg               0.0000
## Avg_superiorfrontal_thickavg                    0.0000
## Avg_superiorparietal_thickavg                   0.0000
## Avg_superiortemporal_thickavg                   0.0000
## Avg_supramarginal_thickavg                      0.0000
## Avg_frontalpole_thickavg                        0.0000
## Avg_temporalpole_thickavg                       0.0000
## Avg_transversetemporal_thickavg                 0.0154
## Avg_insula_thickavg                             0.0000
## Avg_bankssts_surfavg                            1.0000
## Avg_caudalanteriorcingulate_surfavg             0.0000
## Avg_caudalmiddlefrontal_surfavg                 1.0000
## Avg_cuneus_surfavg                              1.0000
## Avg_entorhinal_surfavg                          1.0000
## Avg_fusiform_surfavg                            0.0000
## Avg_inferiorparietal_surfavg                    1.0000
## Avg_inferiortemporal_surfavg                    1.0000
## Avg_isthmuscingulate_surfavg                    1.0000
## Avg_lateraloccipital_surfavg                    1.0000
## Avg_lateralorbitofrontal_surfavg                1.0000
## Avg_lingual_surfavg                             1.0000
## Avg_medialorbitofrontal_surfavg                 1.0000
## Avg_middletemporal_surfavg                      1.0000
## Avg_parahippocampal_surfavg                     0.3080
## Avg_paracentral_surfavg                         1.0000
## Avg_parsopercularis_surfavg                     0.4158
## Avg_parsorbitalis_surfavg                       1.0000
## Avg_parstriangularis_surfavg                    1.0000
## Avg_pericalcarine_surfavg                       1.0000
## Avg_postcentral_surfavg                         1.0000
## Avg_posteriorcingulate_surfavg                  0.8239
## Avg_precentral_surfavg                          1.0000
## Avg_precuneus_surfavg                           1.0000
## Avg_rostralanteriorcingulate_surfavg            0.6083
## Avg_rostralmiddlefrontal_surfavg                0.4158
## Avg_superiorfrontal_surfavg                     1.0000
## Avg_superiorparietal_surfavg                    1.0000
## Avg_superiortemporal_surfavg                    1.0000
## Avg_supramarginal_surfavg                       1.0000
## Avg_frontalpole_surfavg                         1.0000
## Avg_temporalpole_surfavg                        1.0000
## Avg_transversetemporal_surfavg                  0.0077
## Avg_insula_surfavg                              1.0000
## ICV                                             1.0000
## Avg_LatVent                                     0.0000
## Avg_thal                                        0.0000
## Avg_caud                                        0.0000
## Avg_put                                         0.0077
## Avg_pal                                         0.0000
## Avg_hippo                                       0.0000
## Avg_amyg                                        0.0000
## Avg_accumb                                      0.0077
##                                       bonf.age.sex.edu.adj.p.value
## Avg_bankssts_thickavg                                       0.0000
## Avg_caudalanteriorcingulate_thickavg                        1.0000
## Avg_caudalmiddlefrontal_thickavg                            0.0000
## Avg_cuneus_thickavg                                         0.2079
## Avg_entorhinal_thickavg                                     0.6622
## Avg_fusiform_thickavg                                       0.0000
## Avg_inferiorparietal_thickavg                               0.0000
## Avg_inferiortemporal_thickavg                               0.0000
## Avg_isthmuscingulate_thickavg                               0.0000
## Avg_lateraloccipital_thickavg                               0.0000
## Avg_lateralorbitofrontal_thickavg                           0.0000
## Avg_lingual_thickavg                                        1.0000
## Avg_medialorbitofrontal_thickavg                            0.0000
## Avg_middletemporal_thickavg                                 0.0000
## Avg_parahippocampal_thickavg                                0.0924
## Avg_paracentral_thickavg                                    0.0000
## Avg_parsopercularis_thickavg                                0.0000
## Avg_parsorbitalis_thickavg                                  0.0000
## Avg_parstriangularis_thickavg                               0.0000
## Avg_pericalcarine_thickavg                                  1.0000
## Avg_postcentral_thickavg                                    0.0000
## Avg_posteriorcingulate_thickavg                             0.0000
## Avg_precentral_thickavg                                     0.0000
## Avg_precuneus_thickavg                                      0.0000
## Avg_rostralanteriorcingulate_thickavg                       0.0000
## Avg_rostralmiddlefrontal_thickavg                           0.0000
## Avg_superiorfrontal_thickavg                                0.0000
## Avg_superiorparietal_thickavg                               0.0000
## Avg_superiortemporal_thickavg                               0.0000
## Avg_supramarginal_thickavg                                  0.0000
## Avg_frontalpole_thickavg                                    0.0000
## Avg_temporalpole_thickavg                                   0.0000
## Avg_transversetemporal_thickavg                             1.0000
## Avg_insula_thickavg                                         0.0000
## Avg_bankssts_surfavg                                        1.0000
## Avg_caudalanteriorcingulate_surfavg                         0.0000
## Avg_caudalmiddlefrontal_surfavg                             1.0000
## Avg_cuneus_surfavg                                          1.0000
## Avg_entorhinal_surfavg                                      1.0000
## Avg_fusiform_surfavg                                        0.0693
## Avg_inferiorparietal_surfavg                                1.0000
## Avg_inferiortemporal_surfavg                                1.0000
## Avg_isthmuscingulate_surfavg                                1.0000
## Avg_lateraloccipital_surfavg                                1.0000
## Avg_lateralorbitofrontal_surfavg                            1.0000
## Avg_lingual_surfavg                                         1.0000
## Avg_medialorbitofrontal_surfavg                             1.0000
## Avg_middletemporal_surfavg                                  1.0000
## Avg_parahippocampal_surfavg                                 1.0000
## Avg_paracentral_surfavg                                     1.0000
## Avg_parsopercularis_surfavg                                 1.0000
## Avg_parsorbitalis_surfavg                                   1.0000
## Avg_parstriangularis_surfavg                                1.0000
## Avg_pericalcarine_surfavg                                   1.0000
## Avg_postcentral_surfavg                                     1.0000
## Avg_posteriorcingulate_surfavg                              1.0000
## Avg_precentral_surfavg                                      1.0000
## Avg_precuneus_surfavg                                       1.0000
## Avg_rostralanteriorcingulate_surfavg                        0.6622
## Avg_rostralmiddlefrontal_surfavg                            1.0000
## Avg_superiorfrontal_surfavg                                 1.0000
## Avg_superiorparietal_surfavg                                1.0000
## Avg_superiortemporal_surfavg                                1.0000
## Avg_supramarginal_surfavg                                   1.0000
## Avg_frontalpole_surfavg                                     1.0000
## Avg_temporalpole_surfavg                                    1.0000
## Avg_transversetemporal_surfavg                              0.1617
## Avg_insula_surfavg                                          1.0000
## ICV                                                         1.0000
## Avg_LatVent                                                 0.0154
## Avg_thal                                                    0.0000
## Avg_caud                                                    1.0000
## Avg_put                                                     1.0000
## Avg_pal                                                     0.0000
## Avg_hippo                                                   0.0000
## Avg_amyg                                                    0.0000
## Avg_accumb                                                  1.0000
# add sex in training
Train_HC_data$Sex <- as.numeric(Train_HC_data$Sex)
Test_HC_data$Sex <- as.numeric(Test_HC_data$Sex)
Test_Patient_dat$Sex <- as.numeric(Test_Patient_dat$Sex)
Test_schizophrenia_dat$Sex <- as.numeric(Test_schizophrenia_dat$Sex)
Test_TRS_dat$Sex <- as.numeric(Test_TRS_dat$Sex)
Test_FESSD_dat$Sex <- as.numeric(Test_FESSD_dat$Sex)


Ridge (Nested CV)


library(caret)
## 필요한 패키지를 로딩중입니다: lattice
library(parallel)

Nfolds <- 10

f <- as.formula(paste("Age~Sex+",paste(colnames(Train_HC_data)[5:81],collapse = "+"),sep=""))
fitcontrol <- trainControl(method = "repeatedcv", number=Nfolds, search = "random")

ridge.fun <- function(N,i){
  Nfolds <- N
  set.seed(i)
  outer_folds <- createFolds(Train_HC_data$Age,k=Nfolds)
  cv.glmnet <- lapply(outer_folds,function(index) {
    tr <- Train_HC_data[-index,]
    ts <- Train_HC_data[index,]
    fit <- train(f,data=tr,method="ridge",metric="MAE",
               trControl=fitcontrol)
    c(postResample(predict(fit,ts),ts$Age),lambda=fit$bestTune$lambda)
    })
}

Ntimes <- 10
RES <- NULL
for(i in 1:Ntimes){
  res <- mclapply(10,i,FUN=ridge.fun,mc.cores = 6)
  RES <- rbind(RES,t(data.frame(res)))
}
 
library(glmnet)
## 필요한 패키지를 로딩중입니다: Matrix
## Loaded glmnet 4.1-6
ridge.fit <- glmnet(as.matrix(Train_HC_data[,c(3,5:81)]),Train_HC_data[,2],
                    alpha=0,lambda = RES[,4][which.min(RES[,3])])




# Train Control

ctrl_pred_age <- predict(ridge.fit, newx = as.matrix(Train_HC_data[,c(3,5:81)]))
train_error_ctrl_rid <- ctrl_pred_age-Train_HC_data$Age

par(mfrow=c(1,1))
plot(ctrl_pred_age,Train_HC_data$Age,
     xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Train control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
       legend=c(paste("r=",round(cor(Train_HC_data$Age,ctrl_pred_age),4)),
                paste("R2=",round(cor(Train_HC_data$Age,ctrl_pred_age)^2,4)),
                paste("MAE=",round(sum(abs(train_error_ctrl_rid)/length(train_error_ctrl_rid)),4)),
                paste("VEcv=",round(1-sum(train_error_ctrl_rid^2/(var(Train_HC_data$Age)*(length(Train_HC_data$Age)-1))),4)*100,"%")) ,bty="n")

# Test Control

ctrl_pred_age <- predict(ridge.fit, newx = as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_ctrl_rid <- ctrl_pred_age-Test_HC_data$Age

plot(ctrl_pred_age,Test_HC_data$Age,
     xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Test control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
       legend=c(paste("r=",round(cor(Test_HC_data$Age,ctrl_pred_age),4)),
                paste("R2=",round(cor(Test_HC_data$Age,ctrl_pred_age)^2,4)),
                paste("MAE=",round(sum(abs(test_error_ctrl_rid)/length(test_error_ctrl_rid)),4)),
                paste("VEcv=",round(1-sum(test_error_ctrl_rid^2/(var(Test_HC_data$Age)*(length(Test_HC_data$Age)-1))),4)*100,"%")) ,bty="n")

# Patient

pat_pred_age <- predict(ridge.fit, newx = as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_pat_rid <- pat_pred_age-Test_Patient_dat$Age

plot(pat_pred_age,Test_Patient_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Patient)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_Patient_dat$Age,pat_pred_age),4)),
                          paste("R2=",round(cor(Test_Patient_dat$Age,pat_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_pat_rid)/length(test_error_pat_rid)),4)),
paste("VEcv=",round(1-sum(test_error_pat_rid^2/(var(Test_Patient_dat$Age)*(length(Test_Patient_dat$Age)-1))),4)*100,"%")),bty="n")

# Patient - Schizo
sch_pred_age <- predict(ridge.fit, newx = as.matrix(Test_schizophrenia_dat[,c(3,5:81)]))
test_error_sch_rid <- sch_pred_age-Test_schizophrenia_dat$Age

plot(sch_pred_age,Test_schizophrenia_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Schizophrenia)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age),4)),
                          paste("R2=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_sch_rid)/length(test_error_sch_rid)),4)),
paste("VEcv=",round(1-sum(test_error_sch_rid^2/(var(Test_schizophrenia_dat$Age)*(length(Test_schizophrenia_dat$Age)-1))),4)*100,"%")),bty="n")

# Patient - TRS
TRS_pred_age <- predict(ridge.fit, newx = as.matrix(Test_TRS_dat[,c(3,5:81)]))
test_error_trs_rid <- TRS_pred_age-Test_TRS_dat$Age

plot(TRS_pred_age,Test_TRS_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (TRS)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_TRS_dat$Age,TRS_pred_age),4)),
                          paste("R2=",round(cor(Test_TRS_dat$Age,TRS_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_trs_rid)/length(test_error_trs_rid)),4)),
paste("VEcv=",round(1-sum(test_error_trs_rid^2/(var(Test_TRS_dat$Age)*(length(Test_TRS_dat$Age)-1))),4)*100,"%")),bty="n")

# Patient - FESSD
SSD_pred_age <- predict(ridge.fit, newx = as.matrix(Test_FESSD_dat[,c(3,5:81)]))
test_error_ssd_rid <- SSD_pred_age-Test_FESSD_dat$Age

plot(SSD_pred_age,Test_FESSD_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (FE-SSD)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_FESSD_dat$Age,SSD_pred_age),4)),
                          paste("R2=",round(cor(Test_FESSD_dat$Age,SSD_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_ssd_rid)/length(test_error_ssd_rid)),4)),
paste("VEcv=",round(1-sum(test_error_ssd_rid^2/(var(Test_FESSD_dat$Age)*(length(Test_FESSD_dat$Age)-1))),4)*100,"%")),bty="n")

Brain PAD correction-plot

#### Brain PAD correction-HC
pred_age <- predict(ridge.fit,newx=as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_rid <- pred_age-Test_HC_data$Age
PAD_rid <- test_error_rid
cor.test(PAD_rid,Test_HC_data$Age)$p.value
## [1] 2.735224e-17
png("1_rigde_corr.png", width = 800, height = 600)
par(mfrow=c(2,2), oma=c(1,3,1,1))

plot(Test_HC_data$Age,PAD_rid,xlim=c(10,80),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test HC)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_rid~Test_HC_data$Age))
legend("topright",legend=paste("r=",round(cor(PAD_rid,Test_HC_data$Age),4),"(p<0.05)",sep=""))
mtext("A", side=3, line=-3, cex=2, adj = -0.0, outer=TRUE)


Model_linear <- lm(PAD_rid~Test_HC_data$Age+Test_HC_data$Sex)
NullModel <- lm(PAD_rid~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ Test_HC_data$Age + Test_HC_data$Sex
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    208 16868                           
## 2    206 11881  2    4987.8 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Test_HC_data$Age2 <- Test_HC_data$Age^2

Model_quad <- lm(PAD_rid~Test_HC_data$Sex+Test_HC_data$Age+Test_HC_data$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    208 16868                           
## 2    205 11818  3    5050.3 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rid ~ Test_HC_data$Age + Test_HC_data$Sex
## Model 2: PAD_rid ~ Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
##   Res.Df   RSS Df Sum of Sq Pr(>Chi)
## 1    206 11881                      
## 2    205 11818  1    62.413   0.2981
cPAD_rid <- lm(PAD_rid~+Test_HC_data$Sex+Test_HC_data$Age)$residuals

plot(Test_HC_data$Age,cPAD_rid,xlim=c(10,75),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test HC)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_rid~Test_HC_data$Age))
abline(h=0,lty=2,col="red")
mtext("B", side=3, line=-3, cex=2, adj = 0.5, outer=TRUE)


#### Brain PAD correction-Patient
pred_age <- predict(ridge.fit,newx=as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_rid <- pred_age-Test_Patient_dat$Age
PAD_rid <- test_error_rid

plot(Test_Patient_dat$Age,PAD_rid,xlim=c(10,80),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test Patient)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_rid~Test_Patient_dat$Age))
legend("topright",legend=paste("r=",round(cor(PAD_rid,Test_Patient_dat$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_rid,Test_Patient_dat$Age)$p.value
## [1] 3.574547e-10
mtext("C", side=3, line=-27, cex=2, adj = -0.0, outer=TRUE)


Model_linear <- lm(PAD_rid~Test_Patient_dat$Age+Test_Patient_dat$Sex)
NullModel <- lm(PAD_rid~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
##   Res.Df   RSS Df Sum of Sq Pr(>Chi)    
## 1    232 24702                          
## 2    230 20756  2      3946  3.2e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Test_Patient_dat$Age2 <- Test_Patient_dat$Age^2

Model_quad <- lm(PAD_rid~Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    232 24702                           
## 2    229 20335  3    4367.5 1.191e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Model_linear,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rid ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
## Model 2: PAD_rid ~ Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
##   Res.Df   RSS Df Sum of Sq Pr(>Chi)  
## 1    230 20756                        
## 2    229 20335  1    421.53  0.02935 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cPAD_rid <- lm(PAD_rid~+Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)$residuals

plot(Test_Patient_dat$Age,cPAD_rid,xlim=c(10,75),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test Patient)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_rid~Test_Patient_dat$Age))
abline(h=0,lty=2,col="red")
mtext("D", side=3, line=-27, cex=2, adj = 0.5, outer=TRUE)

dev.off()
## png 
##   2

Brain PAD correction -table

test_dat <- rbind(Test_HC_data,Test_Patient_dat)
test_dat$Dx <- c(rep(0,nrow(Test_HC_data)),rep(1,nrow(Test_Patient_dat)))
test_dat$Dx <- factor(test_dat$Dx)

pred_age <- predict(ridge.fit,newx=as.matrix(test_dat[,c(3,5:81)]))
test_error_rid <- pred_age-test_dat$Age
PAD_rid <- test_error_rid
mean(PAD_rid) ; sd(PAD_rid) ; median(PAD_rid)
## [1] 1.513638
## [1] 10.23045
## [1] 2.005724
Model_linear <- lm(PAD_rid~test_dat$Age+test_dat$Sex)
NullModel <- lm(PAD_rid~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ test_dat$Age + test_dat$Sex
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    441 46156                           
## 2    439 36892  2    9264.3 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
test_dat$Age2 <- test_dat$Age^2

Model_quad <- lm(PAD_rid~test_dat$Sex+test_dat$Age+test_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rid ~ 1
## Model 2: PAD_rid ~ test_dat$Sex + test_dat$Age + test_dat$Age2
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    441 46156                           
## 2    438 36157  3    9999.4 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Model_linear,Model_quad,test="Chisq") # Reject H0 
## Analysis of Variance Table
## 
## Model 1: PAD_rid ~ test_dat$Age + test_dat$Sex
## Model 2: PAD_rid ~ test_dat$Sex + test_dat$Age + test_dat$Age2
##   Res.Df   RSS Df Sum of Sq Pr(>Chi)   
## 1    439 36892                         
## 2    438 36157  1     735.1 0.002844 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cPAD_rid <- lm(PAD_rid~+test_dat$Sex+test_dat$Age+test_dat$Age2)$residuals
mean(cPAD_rid) ; sd(cPAD_rid) ; median(cPAD_rid)
## [1] -4.443404e-16
## [1] 9.054699
## [1] 0.02024823
library(plotrix)
dat <- data.frame(SubjID=test_dat$SubjID,Brain_PAD = c(PAD_rid),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))

dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"   
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

library(dplyr)
palette <- c("FE-SSD"="darkgreen","TRS"="red","Control"="grey","Schizophrenia"="blue","Patient"="orange")

dat %>% dplyr::select(Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))
## Warning: `summarise_each_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `across()` instead.
## ℹ The deprecated feature was likely used in the dplyr package.
##   Please report the issue at <https://github.com/tidyverse/dplyr/issues>.
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## # A tibble: 2 × 4
##   Status   mean    sd std.error
##   <chr>   <dbl> <dbl>     <dbl>
## 1 Control -1.89  9.01     0.623
## 2 Patient  4.56 10.3      0.676
dat %>% dplyr::select(Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_sch     mean    sd std.error
##   <chr>         <dbl> <dbl>     <dbl>
## 1 Control       -1.89  9.01     0.623
## 2 Schizophrenia  4.58 10.5      0.750
## 3 <NA>           4.47  9.42     1.55
dat %>% dplyr::select(Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_trs  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -1.89  9.01     0.623
## 2 TRS         6.97  7.54     1.11 
## 3 <NA>        3.97 10.8      0.792
dat %>% dplyr::select(Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_ssd  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -1.89  9.01     0.623
## 2 FE-SSD      5.19 10.7      1.26 
## 3 <NA>        4.29 10.2      0.801
library(ggplot2)
ggplot(dat[!is.na(dat$status_sch),],aes(x = Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_ssd),],aes(x = Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_trs),],aes(x = Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

summary(lm(Brain_PAD~Status+test_dat$Sex+test_dat$Age+test_dat$Age2,data=dat))
## 
## Call:
## lm(formula = Brain_PAD ~ Status + test_dat$Sex + test_dat$Age + 
##     test_dat$Age2, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.862  -5.379   0.325   5.900  23.850 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.732033   4.157596   0.657   0.5115    
## StatusPatient  6.061460   0.831821   7.287 1.49e-12 ***
## test_dat$Sex  -0.968674   0.847673  -1.143   0.2538    
## test_dat$Age   0.185190   0.213239   0.868   0.3856    
## test_dat$Age2 -0.006135   0.002601  -2.359   0.0188 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.589 on 437 degrees of freedom
## Multiple R-squared:  0.3015, Adjusted R-squared:  0.2951 
## F-statistic: 47.16 on 4 and 437 DF,  p-value: < 2.2e-16
summary(lm(Brain_PAD~status_sch+test_dat$Sex+test_dat$Age+test_dat$Age2,data=dat))
## 
## Call:
## lm(formula = Brain_PAD ~ status_sch + test_dat$Sex + test_dat$Age + 
##     test_dat$Age2, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.265  -5.243   0.229   5.836  23.194 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              3.476826   4.392507   0.792   0.4291    
## status_schSchizophrenia  6.547650   0.892309   7.338 1.22e-12 ***
## test_dat$Sex            -0.806285   0.900654  -0.895   0.3712    
## test_dat$Age             0.127325   0.226187   0.563   0.5738    
## test_dat$Age2           -0.005419   0.002730  -1.985   0.0479 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.649 on 400 degrees of freedom
##   (결측으로 인하여 37개의 관측치가 삭제되었습니다.)
## Multiple R-squared:  0.2978, Adjusted R-squared:  0.2908 
## F-statistic: 42.41 on 4 and 400 DF,  p-value: < 2.2e-16
dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]
##     SubjID Brain_PAD  Status status_ssd status_trs    status_sch
## 415    T19 24.279541 Patient     FE-SSD        TRS Schizophrenia
## 419    T24  3.990798 Patient     FE-SSD        TRS Schizophrenia
dat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA

table(dat_for_regression$status_trs)
## 
## Control     TRS 
##     209      46
table(dat_for_regression$status_ssd)
## 
## Control  FE-SSD 
##     209      72
table(dat_for_regression$status_sch)
## 
##       Control Schizophrenia 
##           209           198
test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])

dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]

summary(lm(Brain_PAD~dat_for_regression$status_ssd+test_dat_for_regression$Sex+
             test_dat_for_regression$Age+
             test_dat_for_regression$Age2,data=dat_for_regression))
## 
## Call:
## lm(formula = Brain_PAD ~ dat_for_regression$status_ssd + test_dat_for_regression$Sex + 
##     test_dat_for_regression$Age + test_dat_for_regression$Age2, 
##     data = dat_for_regression)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.0328  -5.2602  -0.3887   5.6345  24.1575 
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                          7.579598   4.626825   1.638 0.102363    
## dat_for_regression$status_ssdFE-SSD  4.295696   1.101383   3.900 0.000117 ***
## dat_for_regression$status_ssdTRS    10.094969   1.366464   7.388  1.3e-12 ***
## test_dat_for_regression$Sex         -0.478314   0.912560  -0.524 0.600539    
## test_dat_for_regression$Age         -0.077507   0.249808  -0.310 0.756560    
## test_dat_for_regression$Age2        -0.003534   0.003025  -1.169 0.243448    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.753 on 321 degrees of freedom
##   (결측으로 인하여 117개의 관측치가 삭제되었습니다.)
## Multiple R-squared:  0.4011, Adjusted R-squared:  0.3917 
## F-statistic: 42.99 on 5 and 321 DF,  p-value: < 2.2e-16
# ANOVA results (Control vs FE-SSD vs TRS)
summary(res_aov <- aov(dat_for_regression$Brain_PAD ~ dat_for_regression$status_ssd))
##                                Df Sum Sq Mean Sq F value   Pr(>F)    
## dat_for_regression$status_ssd   2   4642  2321.2   27.28 1.13e-11 ***
## Residuals                     324  27572    85.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.
pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "fdr")$p.value
##             Control    FE-SSD
## FE-SSD 6.466828e-08        NA
## TRS    2.793844e-08 0.3059112
pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value
##             Control    FE-SSD
## FE-SSD 1.293366e-07        NA
## TRS    2.793844e-08 0.9177336
# t.test results (Control vs Schizo)
t.test(dat$Brain_PAD~dat$status_sch)
## 
##  Welch Two Sample t-test
## 
## data:  dat$Brain_PAD by dat$status_sch
## t = -6.6336, df = 385.01, p-value = 1.112e-10
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -8.385513 -4.551188
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -1.887131                    4.581220
t.test(dat$Brain_PAD~dat$status_sch, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  dat$Brain_PAD by dat$status_sch
## t = -6.6663, df = 403, p-value = 8.657e-11
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -8.375855 -4.560847
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -1.887131                    4.581220
#### corrected_Brain_PAD ####
dat <- data.frame(SubjID=test_dat$SubjID,corrected_Brain_PAD = c(cPAD_rid),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))

dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"   
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

dat %>% dplyr::select(corrected_Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 2 × 4
##   Status   mean    sd std.error
##   <chr>   <dbl> <dbl>     <dbl>
## 1 Control -3.09  7.57     0.524
## 2 Patient  2.77  9.39     0.615
dat %>% dplyr::select(corrected_Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_sch      mean    sd std.error
##   <chr>          <dbl> <dbl>     <dbl>
## 1 Control       -3.09   7.57     0.524
## 2 Schizophrenia  3.11   9.65     0.690
## 3 <NA>           0.993  7.69     1.26
dat %>% dplyr::select(corrected_Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_trs  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -3.09  7.57     0.524
## 2 TRS         5.77  7.32     1.08 
## 3 <NA>        2.04  9.70     0.710
dat %>% dplyr::select(corrected_Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_ssd  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -3.09  7.57     0.524
## 2 FE-SSD      1.85  8.66     1.02 
## 3 <NA>        3.19  9.69     0.764
library(ggplot2)
ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette[c(3,5)]) +
  coord_cartesian(xlim =c(-30,30)) + theme_bw()

A <- ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette[c(3,4)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

B <- ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette[c(3,1)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

C <- ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette[c(3,2)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

png("S5-2(ridge).png", width = 1400, height = 300)
ggarrange(A,B,C,
          labels = c("A","B","C"), 
          ncol = 3, nrow = 1,
          font.label = list(size=20))
dev.off()
## png 
##   2
# ANOVA results (Control vs FE-SSD vs TRS)
dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]
##     SubjID corrected_Brain_PAD  Status status_ssd status_trs    status_sch
## 981    T19         19.94609021 Patient     FE-SSD        TRS Schizophrenia
## 985    T24         -0.03676607 Patient     FE-SSD        TRS Schizophrenia
dat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA

table(dat_for_regression$status_trs)
## 
## Control     TRS 
##     209      46
table(dat_for_regression$status_ssd)
## 
## Control  FE-SSD 
##     209      72
test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])

dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]

summary(res_aov <- aov(dat_for_regression$corrected_Brain_PAD~
                         dat_for_regression$status_ssd))
##                                Df Sum Sq Mean Sq F value   Pr(>F)    
## dat_for_regression$status_ssd   2   3590  1794.8   29.58 1.59e-12 ***
## Residuals                     324  19661    60.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.
pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
               dat_for_regression$status_ssd, p.adj = "fdr")$p.value
##             Control      FE-SSD
## FE-SSD 7.577113e-06          NA
## TRS    4.792668e-11 0.007980272
pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
               dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value
##             Control     FE-SSD
## FE-SSD 1.515423e-05         NA
## TRS    4.792668e-11 0.02394082
# t.test results (Control vs Schizo)
t.test(dat$corrected_Brain_PAD~dat$status_sch)
## 
##  Welch Two Sample t-test
## 
## data:  dat$corrected_Brain_PAD by dat$status_sch
## t = -7.1625, df = 369.56, p-value = 4.317e-12
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -7.904795 -4.499344
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -3.092237                    3.109832
t.test(dat$corrected_Brain_PAD~dat$status_sch, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  dat$corrected_Brain_PAD by dat$status_sch
## t = -7.2177, df = 403, p-value = 2.65e-12
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -7.891324 -4.512815
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -3.092237                    3.109832
summary(res_aov <- aov(dat$corrected_Brain_PAD~dat$status_sch))
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## dat$status_sch   1   3891    3891   52.09 2.65e-12 ***
## Residuals      403  30098      75                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 37개의 관측치가 삭제되었습니다.


SVR[RBF] (Nested CV)


library(kernlab)
## 
## 다음의 패키지를 부착합니다: 'kernlab'
## The following object is masked from 'package:ggplot2':
## 
##     alpha
start.time <- Sys.time()

Nfolds <- 10

svm.fun <- function(N,i){
  Nfolds <- N
  set.seed(i)
  outer_folds <- createFolds(Train_HC_data$Age,k=Nfolds)
  cv.svm <- lapply(outer_folds,function(index) {
    tr <- Train_HC_data[-index,]
    ts <- Train_HC_data[index,]
    svm.fit <- ksvm(as.matrix(tr[, c(3,5:81)]),
               tr[,2],
               kpar="automatic",
               kernel="rbfdot",cross=Nfolds) 
    
    c(svm.fit@cross,svm.fit@kernelf@kpar[["sigma"]])
    })
}

Ntimes <- 10
RES <- NULL
for(i in 1:Ntimes){
  res <- mclapply(10,i,FUN=svm.fun,mc.cores = 6)
  RES <- rbind(RES,t(data.frame(res)))
}
 
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken
## Time difference of 29.7719 secs
svm.fit.77 <- ksvm(as.matrix(Train_HC_data[,c(3,5:81)]),Train_HC_data[,2],
                  kpar=list(sigma=RES[,2][which.min(RES[,1])]),
                  kernel="rbfdot")
svm.fit.77
## Support Vector Machine object of class "ksvm" 
## 
## SV type: eps-svr  (regression) 
##  parameter : epsilon = 0.1  cost C = 1 
## 
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  0.0077714944308448 
## 
## Number of Support Vectors : 449 
## 
## Objective Function Value : -144.8082 
## Training error : 0.141236
# Train Control

ctrl_pred_age <- predict(svm.fit.77 ,as.matrix(Train_HC_data[,c(3,5:81)]))
train_error_ctrl_svm <- ctrl_pred_age-Train_HC_data$Age

par(mfrow=c(1,1))
plot(ctrl_pred_age,Train_HC_data$Age,
     xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Train control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
       legend=c(paste("r=",round(cor(Train_HC_data$Age,ctrl_pred_age),4)),
                paste("R2=",round(cor(Train_HC_data$Age,ctrl_pred_age)^2,4)),
                paste("MAE=",round(sum(abs(train_error_ctrl_svm)/length(train_error_ctrl_svm)),4)),
                paste("VEcv=",round(1-sum(train_error_ctrl_svm^2/(var(Train_HC_data$Age)*(length(Train_HC_data$Age)-1))),4)*100,"%")) ,bty="n")

# Test Control

ctrl_pred_age <- predict(svm.fit.77 , as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_ctrl_svm <- ctrl_pred_age-Test_HC_data$Age

plot(ctrl_pred_age,Test_HC_data$Age,
     xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Test control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
       legend=c(paste("r=",round(cor(Test_HC_data$Age,ctrl_pred_age),4)),
                paste("R2=",round(cor(Test_HC_data$Age,ctrl_pred_age)^2,4)),
                paste("MAE=",round(sum(abs(test_error_ctrl_svm)/length(test_error_ctrl_svm)),4)),
                paste("VEcv=",round(1-sum(test_error_ctrl_svm^2/(var(Test_HC_data$Age)*(length(Test_HC_data$Age)-1))),4)*100,"%")) ,bty="n")

# Patient
pat_pred_age <- predict(svm.fit.77 ,as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_pat_svm <- pat_pred_age-Test_Patient_dat$Age

plot(pat_pred_age,Test_Patient_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Patient)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_Patient_dat$Age,pat_pred_age),4)),paste("R2=",round(cor(Test_Patient_dat$Age,pat_pred_age)^2,4)),                   paste("MAE=",round(sum(abs(test_error_pat_svm)/length(test_error_pat_svm)),4)),
                          paste("VEcv=",round(1-sum(test_error_pat_svm^2/(var(Test_Patient_dat$Age)*(length(Test_Patient_dat$Age)-1))),4)*100,"%")),bty="n")

# Patient - schizo
sch_pred_age <- predict(svm.fit.77, as.matrix(Test_schizophrenia_dat[,c(3,5:81)]))
test_error_sch_svm <- sch_pred_age-Test_schizophrenia_dat$Age

plot(sch_pred_age,Test_schizophrenia_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Schizophrenia)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age),4)),
                          paste("R2=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_sch_svm)/length(test_error_sch_svm)),4)),
paste("VEcv=",round(1-sum(test_error_sch_svm^2/(var(Test_schizophrenia_dat$Age)*(length(Test_schizophrenia_dat$Age)-1))),4)*100,"%")),bty="n")

# Patient - TRS
TRS_pred_age <- predict(svm.fit.77, as.matrix(Test_TRS_dat[,c(3,5:81)]))
test_error_trs_svm <- TRS_pred_age-Test_TRS_dat$Age

plot(TRS_pred_age,Test_TRS_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (TRS)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_TRS_dat$Age,TRS_pred_age),4)),
                          paste("R2=",round(cor(Test_TRS_dat$Age,TRS_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_trs_svm)/length(test_error_trs_svm)),4)),
paste("VEcv=",round(1-sum(test_error_trs_svm^2/(var(Test_TRS_dat$Age)*(length(Test_TRS_dat$Age)-1))),4)*100,"%")),bty="n")

# Patient - FESSD
SSD_pred_age <- predict(svm.fit.77, as.matrix(Test_FESSD_dat[,c(3,5:81)]))
test_error_ssd_svm <- SSD_pred_age-Test_FESSD_dat$Age

plot(SSD_pred_age,Test_FESSD_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (FE-SSD)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_FESSD_dat$Age,SSD_pred_age),4)),
                          paste("R2=",round(cor(Test_FESSD_dat$Age,SSD_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_ssd_svm)/length(test_error_ssd_svm)),4)),
paste("VEcv=",round(1-sum(test_error_ssd_svm^2/(var(Test_FESSD_dat$Age)*(length(Test_FESSD_dat$Age)-1))),4)*100,"%")),bty="n")

Brain PAD correction-plot

#### Brain PAD correction-HC
pred_age <- predict(svm.fit.77, as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_svm <- pred_age-Test_HC_data$Age
PAD_svm <- test_error_svm

png("2_SVR_corr.png", width = 800, height = 600)
par(mfrow=c(2,2), oma=c(1,3,1,1))

plot(Test_HC_data$Age,PAD_svm,xlim=c(10,80),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test HC)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_svm~Test_HC_data$Age))
legend("topright",legend=paste("r=",round(cor(PAD_svm,Test_HC_data$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_svm,Test_HC_data$Age)$p.value
## [1] 3.178454e-25
mtext("A", side=3, line=-3, cex=2, adj = -0.0, outer=TRUE)


Model_linear <- lm(PAD_svm~Test_HC_data$Age+Test_HC_data$Sex)
NullModel <- lm(PAD_svm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ Test_HC_data$Age + Test_HC_data$Sex
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    208 14952                           
## 2    206  8783  2    6169.3 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Test_HC_data$Age2 <- Test_HC_data$Age^2
Model_quad <- lm(PAD_svm~+Test_HC_data$Sex+Test_HC_data$Age+Test_HC_data$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ +Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
##   Res.Df     RSS Df Sum of Sq  Pr(>Chi)    
## 1    208 14952.3                           
## 2    205  8743.2  3    6209.1 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_svm ~ Test_HC_data$Age + Test_HC_data$Sex
## Model 2: PAD_svm ~ +Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
##   Res.Df    RSS Df Sum of Sq Pr(>Chi)
## 1    206 8783.0                      
## 2    205 8743.2  1    39.742   0.3344
cPAD_svm <- lm(PAD_svm~+Test_HC_data$Sex+Test_HC_data$Age)$residuals


plot(Test_HC_data$Age,cPAD_svm,xlim=c(10,75),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test HC)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_svm~Test_HC_data$Age))
abline(h=0,lty=2,col="red")
mtext("B", side=3, line=-3, cex=2, adj = 0.5, outer=TRUE)


#### Brain PAD correction-Patient
pred_age <- predict(svm.fit.77, as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_svm <- pred_age-Test_Patient_dat$Age
PAD_svm <- test_error_svm


plot(Test_Patient_dat$Age,PAD_svm,xlim=c(10,80),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test Patient)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_svm~Test_Patient_dat$Age))
legend("topright",legend=paste("r=",round(cor(PAD_svm,Test_Patient_dat$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_svm,Test_Patient_dat$Age)$p.value
## [1] 4.532322e-23
mtext("C", side=3, line=-27, cex=2, adj = -0.0, outer=TRUE)

Model_linear <- lm(PAD_svm~Test_Patient_dat$Age+Test_Patient_dat$Sex)
NullModel <- lm(PAD_svm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    232 23999                           
## 2    230 15552  2    8447.1 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Test_Patient_dat$Age2 <- Test_Patient_dat$Age^2
Model_quad <- lm(PAD_svm~+Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ +Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    232 23999                           
## 2    229 14610  3      9389 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Model_linear,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_svm ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
## Model 2: PAD_svm ~ +Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    230 15552                           
## 2    229 14610  1    941.93 0.0001218 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cPAD_svm <- lm(PAD_svm~+Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)$residuals

plot(Test_Patient_dat$Age,cPAD_svm,xlim=c(10,75),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test Patienit)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_svm~Test_Patient_dat$Age+Test_Patient_dat$Age2))
## Warning in abline(lm(cPAD_svm ~ Test_Patient_dat$Age + Test_Patient_dat$Age2)):
## only using the first two of 3 regression coefficients
abline(h=0,lty=2,col="red")
mtext("D", side=3, line=-27, cex=2, adj = 0.5, outer=TRUE)

dev.off()
## png 
##   2

Brain PAD correction -table

test_dat <- rbind(Test_HC_data,Test_Patient_dat)
test_dat$Dx <- c(rep(0,nrow(Test_HC_data)),rep(1,nrow(Test_Patient_dat)))
test_dat$Dx <- factor(test_dat$Dx)

pred_age <- predict(svm.fit.77 , as.matrix(test_dat[,c(3,5:81)]))
test_error_svm <- pred_age-test_dat$Age
PAD_svm <- test_error_svm
mean(PAD_svm) ; sd(PAD_svm) ; median(PAD_svm)
## [1] 0.07811399
## [1] 9.6726
## [1] 0.6968941
Model_linear <- lm(PAD_svm~test_dat$Age+test_dat$Sex)
NullModel <- lm(PAD_svm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ test_dat$Age + test_dat$Sex
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    441 41260                           
## 2    439 26917  2     14342 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
test_dat$Age2 <- test_dat$Age^2
Model_quad <- lm(PAD_svm~+test_dat$Sex+test_dat$Age+test_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_svm ~ 1
## Model 2: PAD_svm ~ +test_dat$Sex + test_dat$Age + test_dat$Age2
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    441 41260                           
## 2    438 26030  3     15230 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_svm ~ test_dat$Age + test_dat$Sex
## Model 2: PAD_svm ~ +test_dat$Sex + test_dat$Age + test_dat$Age2
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    439 26917                           
## 2    438 26030  1    887.43 0.0001114 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cPAD_svm <- lm(PAD_svm~+test_dat$Sex+test_dat$Age+test_dat$Age2)$residuals
mean(cPAD_svm) ; sd(cPAD_svm) ; median(cPAD_svm)
## [1] 2.730541e-16
## [1] 7.682741
## [1] -0.2859405
library(plotrix)
dat <- data.frame(SubjID=test_dat$SubjID,Brain_PAD = c(PAD_svm),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))

dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"   
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

library(dplyr)
palette <- c("FE-SSD"="darkgreen","TRS"="red","Control"="grey","Schizophrenia"="blue","Patient"="orange")

dat %>% dplyr::select(Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 2 × 4
##   Status   mean    sd std.error
##   <chr>   <dbl> <dbl>     <dbl>
## 1 Control -2.33  8.48     0.586
## 2 Patient  2.24 10.2      0.666
dat %>% dplyr::select(Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_sch     mean    sd std.error
##   <chr>         <dbl> <dbl>     <dbl>
## 1 Control       -2.33  8.48     0.586
## 2 Schizophrenia  1.76 10.2      0.728
## 3 <NA>           4.79  9.83     1.62
dat %>% dplyr::select(Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_trs  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -2.33  8.48     0.586
## 2 TRS         3.06  8.11     1.20 
## 3 <NA>        2.04 10.6      0.777
dat %>% dplyr::select(Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_ssd  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -2.33  8.48     0.586
## 2 FE-SSD      4.53 10.8      1.27 
## 3 <NA>        1.22  9.75     0.769
ggplot(dat,aes(x = Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette) +
  coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_sch),],aes(x = Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_ssd),],aes(x = Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_trs),],aes(x = Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

summary(lm(Brain_PAD~Status+test_dat$Sex+test_dat$Age+test_dat$Age2,data=dat))
## 
## Call:
## lm(formula = Brain_PAD ~ Status + test_dat$Sex + test_dat$Age + 
##     test_dat$Age2, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.9510  -4.2134  -0.1874   5.0155  26.1833 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.0002983  3.5769235   0.000 0.999933    
## StatusPatient  4.5084597  0.7156445   6.300 7.28e-10 ***
## test_dat$Sex   1.6882406  0.7292822   2.315 0.021079 *  
## test_dat$Age   0.1881796  0.1834570   1.026 0.305581    
## test_dat$Age2 -0.0074606  0.0022376  -3.334 0.000928 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.39 on 437 degrees of freedom
## Multiple R-squared:  0.4216, Adjusted R-squared:  0.4164 
## F-statistic: 79.65 on 4 and 437 DF,  p-value: < 2.2e-16
summary(lm(Brain_PAD~status_sch+test_dat$Sex+test_dat$Age+test_dat$Age2,data=dat))
## 
## Call:
## lm(formula = Brain_PAD ~ status_sch + test_dat$Sex + test_dat$Age + 
##     test_dat$Age2, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.4697  -4.1100  -0.2096   4.9655  25.9779 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              0.556947   3.737911   0.149  0.88163    
## status_schSchizophrenia  4.728114   0.759332   6.227 1.21e-09 ***
## test_dat$Sex             1.738308   0.766434   2.268  0.02386 *  
## test_dat$Age             0.148296   0.192479   0.770  0.44149    
## test_dat$Age2           -0.006932   0.002323  -2.984  0.00302 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.36 on 400 degrees of freedom
##   (결측으로 인하여 37개의 관측치가 삭제되었습니다.)
## Multiple R-squared:  0.4126, Adjusted R-squared:  0.4067 
## F-statistic: 70.24 on 4 and 400 DF,  p-value: < 2.2e-16
dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]
##     SubjID Brain_PAD  Status status_ssd status_trs    status_sch
## 415    T19 22.432483 Patient     FE-SSD        TRS Schizophrenia
## 419    T24  4.575506 Patient     FE-SSD        TRS Schizophrenia
dat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA

table(dat_for_regression$status_trs)
## 
## Control     TRS 
##     209      46
table(dat_for_regression$status_ssd)
## 
## Control  FE-SSD 
##     209      72
test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])

dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]

summary(lm(Brain_PAD~dat_for_regression$status_ssd+test_dat_for_regression$Sex+
             test_dat_for_regression$Age+
             test_dat_for_regression$Age2,data=dat_for_regression))
## 
## Call:
## lm(formula = Brain_PAD ~ dat_for_regression$status_ssd + test_dat_for_regression$Sex + 
##     test_dat_for_regression$Age + test_dat_for_regression$Age2, 
##     data = dat_for_regression)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.1642  -4.2029  -0.3043   4.9129  20.1372 
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                          5.052667   4.128789   1.224 0.221939    
## dat_for_regression$status_ssdFE-SSD  3.628643   0.982829   3.692 0.000261 ***
## dat_for_regression$status_ssdTRS     7.591653   1.219377   6.226  1.5e-09 ***
## test_dat_for_regression$Sex          1.954197   0.814331   2.400 0.016975 *  
## test_dat_for_regression$Age         -0.110802   0.222918  -0.497 0.619492    
## test_dat_for_regression$Age2        -0.003924   0.002699  -1.454 0.146950    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.918 on 321 degrees of freedom
##   (결측으로 인하여 117개의 관측치가 삭제되었습니다.)
## Multiple R-squared:  0.4736, Adjusted R-squared:  0.4654 
## F-statistic: 57.77 on 5 and 321 DF,  p-value: < 2.2e-16
# ANOVA results (Control vs FE-SSD vs TRS)
summary(res_aov <- aov(dat_for_regression$Brain_PAD ~ dat_for_regression$status_ssd))
##                                Df Sum Sq Mean Sq F value  Pr(>F)    
## dat_for_regression$status_ssd   2   3043  1521.5   18.85 1.8e-08 ***
## Residuals                     324  26146    80.7                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.
pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "fdr")$p.value
##             Control    FE-SSD
## FE-SSD 1.453155e-07        NA
## TRS    4.012511e-04 0.3868526
pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value
##             Control FE-SSD
## FE-SSD 1.453155e-07     NA
## TRS    8.025021e-04      1
# t.test results (Control vs Schizo)
t.test(dat$Brain_PAD~dat$status_sch)
## 
##  Welch Two Sample t-test
## 
## data:  dat$Brain_PAD by dat$status_sch
## t = -4.383, df = 380.21, p-value = 1.516e-05
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -5.933928 -2.258724
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -2.334937                    1.761390
t.test(dat$Brain_PAD~dat$status_sch, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  dat$Brain_PAD by dat$status_sch
## t = -4.4087, df = 403, p-value = 1.335e-05
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -5.922893 -2.269759
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -2.334937                    1.761390
#### corrected_Brain_PAD ####

dat <- data.frame(SubjID=test_dat$SubjID,corrected_Brain_PAD = c(cPAD_svm),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))

dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"   
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

dat %>% dplyr::select(corrected_Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 2 × 4
##   Status   mean    sd std.error
##   <chr>   <dbl> <dbl>     <dbl>
## 1 Control -2.30  6.59     0.456
## 2 Patient  2.06  8.02     0.525
dat %>% dplyr::select(corrected_Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_sch     mean    sd std.error
##   <chr>         <dbl> <dbl>     <dbl>
## 1 Control       -2.30  6.59     0.456
## 2 Schizophrenia  2.19  8.08     0.577
## 3 <NA>           1.40  7.74     1.27
dat %>% dplyr::select(corrected_Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_trs  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -2.30  6.59     0.456
## 2 TRS         4.34  6.32     0.932
## 3 <NA>        1.50  8.30     0.607
dat %>% dplyr::select(corrected_Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_ssd  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -2.30  6.59     0.456
## 2 FE-SSD      1.67  8.18     0.964
## 3 <NA>        2.24  7.96     0.627
ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette) +
  coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

library(ggplot2)
ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette[c(3,5)]) +
  coord_cartesian(xlim =c(-30,30)) + theme_bw()

A <- ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette[c(3,4)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

B <- ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette[c(3,1)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

C <- ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette[c(3,2)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

png("p2(SVR).png", width = 1400, height = 300)
ggarrange(A,B,C,
          labels = c("A","B","C"), 
          ncol = 3, nrow = 1,
          font.label = list(size=20))
dev.off()
## png 
##   2
dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]
##     SubjID corrected_Brain_PAD  Status status_ssd status_trs    status_sch
## 415    T19          16.7968216 Patient     FE-SSD        TRS Schizophrenia
## 419    T24          -0.3460116 Patient     FE-SSD        TRS Schizophrenia
dat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA

table(dat_for_regression$status_trs)
## 
## Control     TRS 
##     209      46
table(dat_for_regression$status_ssd)
## 
## Control  FE-SSD 
##     209      72
test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])

dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]

summary(res_aov <- aov(dat_for_regression$corrected_Brain_PAD~
                         dat_for_regression$status_ssd))
##                                Df Sum Sq Mean Sq F value   Pr(>F)    
## dat_for_regression$status_ssd   2   2094  1046.9   21.77 1.34e-09 ***
## Residuals                     324  15577    48.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.
pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
               dat_for_regression$status_ssd, p.adj = "fdr")$p.value
##             Control     FE-SSD
## FE-SSD 5.380789e-05         NA
## TRS    3.071392e-08 0.04227491
pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
               dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value
##             Control    FE-SSD
## FE-SSD 1.076158e-04        NA
## TRS    3.071392e-08 0.1268247
# t.test results (Control vs Schizo)
t.test(dat$corrected_Brain_PAD~dat$status_sch)
## 
##  Welch Two Sample t-test
## 
## data:  dat$corrected_Brain_PAD by dat$status_sch
## t = -6.1046, df = 376.68, p-value = 2.565e-09
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -5.934865 -3.043082
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -2.299979                    2.188995
t.test(dat$corrected_Brain_PAD~dat$status_sch, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  dat$corrected_Brain_PAD by dat$status_sch
## t = -6.1443, df = 403, p-value = 1.93e-09
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -5.925225 -3.052722
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -2.299979                    2.188995
summary(res_aov <- aov(dat$corrected_Brain_PAD~dat$status_sch))
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## dat$status_sch   1   2038    2038   37.75 1.93e-09 ***
## Residuals      403  21757      54                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 37개의 관측치가 삭제되었습니다.


RVR[RBF] (Nested CV)


start.time <- Sys.time()

Nfolds <- 10

rvm.fun <- function(N,i){
  Nfolds <- N
  set.seed(i)
  outer_folds <- createFolds(Train_HC_data$Age,k=Nfolds)
  cv.rvm <- lapply(outer_folds,function(index) {
    tr <- Train_HC_data[-index,]
    ts <- Train_HC_data[index,]
    rvm.fit <- rvm(as.matrix(tr[, c(3,5:81)]),
               tr[,2],type="regression",
               kpar="automatic",
               kernel="rbfdot",cross=Nfolds)
    
    c(rvm.fit@cross,rvm.fit@kernelf@kpar[["sigma"]])
    })
}

Ntimes <- 10
RES <- NULL
for(i in 1:Ntimes){
  res <- mclapply(10,i,FUN=rvm.fun,mc.cores = 6)
  RES <- rbind(RES,t(data.frame(res)))
}
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel 
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken
## Time difference of 29.99672 mins
rvm.fit.77 <- rvm(as.matrix(Train_HC_data[,c(3,5:81)]),Train_HC_data[,2],
                  kpar=list(sigma=RES[,2][which.min(RES[,1])]),
                  type="regression",kernel="rbfdot")
rvm.fit.77
## Relevance Vector Machine object of class "rvm" 
## Problem type: regression 
##  
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  2.60538261599689e-10 
## 
## Number of Relevance Vectors : 33 
## Variance :  105.0113
## Training error : 99.079268256
# Train Control
ctrl_pred_age <- predict(rvm.fit.77 ,as.matrix(Train_HC_data[,c(3,5:81)]))
train_error_ctrl_rvm <- ctrl_pred_age-Train_HC_data$Age

par(mfrow=c(1,1))
plot(ctrl_pred_age,Train_HC_data$Age,
     xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Train control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
       legend=c(paste("r=",round(cor(Train_HC_data$Age,ctrl_pred_age),4)),
                paste("R2=",round(cor(Train_HC_data$Age,ctrl_pred_age)^2,4)),
                paste("MAE=",round(sum(abs(train_error_ctrl_rid)/length(train_error_ctrl_rid)),4)),
                paste("VEcv=",round(1-sum(train_error_ctrl_rid^2/(var(Train_HC_data$Age)*(length(Train_HC_data$Age)-1))),4)*100,"%")) ,bty="n")

# Test Control

ctrl_pred_age <- predict(rvm.fit.77 , as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_ctrl_rvm <- ctrl_pred_age-Test_HC_data$Age

plot(ctrl_pred_age,Test_HC_data$Age,
     xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Test control data)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",
       legend=c(paste("r=",round(cor(Test_HC_data$Age,ctrl_pred_age),4)),
                paste("R2=",round(cor(Test_HC_data$Age,ctrl_pred_age)^2,4)),
                paste("MAE=",round(sum(abs(test_error_ctrl_rvm)/length(test_error_ctrl_rvm)),4)),
                paste("VEcv=",round(1-sum(test_error_ctrl_rvm^2/(var(Test_HC_data$Age)*(length(Test_HC_data$Age)-1))),4)*100,"%")) ,bty="n")

# Patient

pat_pred_age <- predict(rvm.fit.77 ,as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_pat_rvm <- pat_pred_age-Test_Patient_dat$Age

plot(pat_pred_age,Test_Patient_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Patient)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_Patient_dat$Age,pat_pred_age),4)),paste("R2=",round(cor(Test_Patient_dat$Age,pat_pred_age)^2,4)),
                          paste("MAE=",round(sum(abs(test_error_pat_rid)/length(test_error_pat_rid)),4)),
                          paste("VEcv=",round(1-sum(test_error_pat_rid^2/(var(Test_Patient_dat$Age)*(length(Test_Patient_dat$Age)-1))),4)*100,"%")),bty="n")

# Patient - schizo
sch_pred_age <- predict(rvm.fit.77, as.matrix(Test_schizophrenia_dat[,c(3,5:81)]))
test_error_sch_rvm <- sch_pred_age-Test_schizophrenia_dat$Age

plot(sch_pred_age,Test_schizophrenia_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (Schizophrenia)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age),4)),
                          paste("R2=",round(cor(Test_schizophrenia_dat$Age,sch_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_sch_rvm)/length(test_error_sch_rvm)),4)),
paste("VEcv=",round(1-sum(test_error_sch_rvm^2/(var(Test_schizophrenia_dat$Age)*(length(Test_schizophrenia_dat$Age)-1))),4)*100,"%")),bty="n")

# Patient - TRS
TRS_pred_age <- predict(rvm.fit.77, as.matrix(Test_TRS_dat[,c(3,5:81)]))
test_error_trs_rvm <- TRS_pred_age-Test_TRS_dat$Age

plot(TRS_pred_age,Test_TRS_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (TRS)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_TRS_dat$Age,TRS_pred_age),4)),
                          paste("R2=",round(cor(Test_TRS_dat$Age,TRS_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_trs_rvm)/length(test_error_trs_rvm)),4)),
paste("VEcv=",round(1-sum(test_error_trs_rvm^2/(var(Test_TRS_dat$Age)*(length(Test_TRS_dat$Age)-1))),4)*100,"%")),bty="n")

# Patient - FESSD
SSD_pred_age <- predict(rvm.fit.77, as.matrix(Test_FESSD_dat[,c(3,5:81)]))
test_error_ssd_rvm <- SSD_pred_age-Test_FESSD_dat$Age

plot(SSD_pred_age,Test_FESSD_dat$Age,xlim=c(10,80),ylim=c(10,80),ylab="Chronological age (FE-SSD)",xlab="Predicted brain age",pch=16)
abline(0,1)
legend("topleft",legend=c(paste("r=",round(cor(Test_FESSD_dat$Age,SSD_pred_age),4)),
                          paste("R2=",round(cor(Test_FESSD_dat$Age,SSD_pred_age)^2,4)),
paste("MAE=",round(sum(abs(test_error_ssd_rvm)/length(test_error_ssd_rvm)),4)),
paste("VEcv=",round(1-sum(test_error_ssd_rvm^2/(var(Test_FESSD_dat$Age)*(length(Test_FESSD_dat$Age)-1))),4)*100,"%")),bty="n")

Brain PAD correction-plot

#### Brain PAD correction-HC
pred_age <- predict(rvm.fit.77 , as.matrix(Test_HC_data[,c(3,5:81)]))
test_error_rvm <- pred_age-Test_HC_data$Age
PAD_rvm <- test_error_rvm

png("3_RVR_corr.png", width = 800, height = 600)
par(mfrow=c(2,2), oma=c(1,3,1,1))

plot(Test_HC_data$Age,PAD_rvm,xlim=c(10,80),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test HC)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_rvm~Test_HC_data$Age))
legend("topright",legend=paste("r=",round(cor(PAD_rvm,Test_HC_data$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_rvm,Test_HC_data$Age)$p.value
## [1] 9.523509e-56
mtext("A", side=3, line=-3, cex=2, adj = -0.0, outer=TRUE)

Model_linear <- lm(PAD_rvm~Test_HC_data$Age+Test_HC_data$Sex)
NullModel <- lm(PAD_rvm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ Test_HC_data$Age + Test_HC_data$Sex
##   Res.Df     RSS Df Sum of Sq  Pr(>Chi)    
## 1    208 30264.1                           
## 2    206  8630.2  2     21634 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Test_HC_data$Age2 <- Test_HC_data$Age^2
Model_quad <- lm(PAD_rvm~Test_HC_data$Sex+Test_HC_data$Age+Test_HC_data$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
##   Res.Df     RSS Df Sum of Sq  Pr(>Chi)    
## 1    208 30264.1                           
## 2    205  8451.4  3     21813 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Model_linear,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rvm ~ Test_HC_data$Age + Test_HC_data$Sex
## Model 2: PAD_rvm ~ Test_HC_data$Sex + Test_HC_data$Age + Test_HC_data$Age2
##   Res.Df    RSS Df Sum of Sq Pr(>Chi)  
## 1    206 8630.2                        
## 2    205 8451.4  1    178.74  0.03732 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# cPAD_rvm <- lm(PAD_rvm~Test_HC_data$Sex+Test_HC_data$Age)$residuals
cPAD_rvm <- lm(PAD_rvm~Test_HC_data$Sex+Test_HC_data$Age+Test_HC_data$Age2)$residuals

plot(Test_HC_data$Age,cPAD_rvm,xlim=c(10,75),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test HC)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_rvm~Test_HC_data$Age))
abline(h=0,lty=2,col="red")
mtext("B", side=3, line=-3, cex=2, adj = 0.5, outer=TRUE)



#### Brain PAD correction-Patient
pred_age <- predict(rvm.fit.77 , as.matrix(Test_Patient_dat[,c(3,5:81)]))
test_error_rvm <- pred_age-Test_Patient_dat$Age
PAD_rvm <- test_error_rvm

plot(Test_Patient_dat$Age,PAD_rvm,xlim=c(10,80),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test Patient)",ylab="Brain PAD uncorrected")
abline(h=0,lty=2,col="red")
abline(lm(PAD_rvm~Test_Patient_dat$Age))
legend("topright",legend=paste("r=",round(cor(PAD_rvm,Test_Patient_dat$Age),4),"(p<0.05)",sep=""))
cor.test(PAD_rvm,Test_Patient_dat$Age)$p.value
## [1] 7.696674e-32
mtext("C", side=3, line=-27, cex=2, adj = -0.0, outer=TRUE)

Model_linear <- lm(PAD_rvm~Test_Patient_dat$Age+Test_Patient_dat$Sex)
NullModel <- lm(PAD_rvm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    232 33880                           
## 2    230 17664  2     16216 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Test_Patient_dat$Age2 <- Test_Patient_dat$Age^2
Model_quad <- lm(PAD_rvm~Test_Patient_dat$Sex+Test_Patient_dat$Age+Test_Patient_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    232 33880                           
## 2    229 17628  3     16252 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rvm ~ Test_Patient_dat$Age + Test_Patient_dat$Sex
## Model 2: PAD_rvm ~ Test_Patient_dat$Sex + Test_Patient_dat$Age + Test_Patient_dat$Age2
##   Res.Df   RSS Df Sum of Sq Pr(>Chi)
## 1    230 17664                      
## 2    229 17628  1    36.347    0.492
cPAD_rvm <- lm(PAD_rvm~Test_Patient_dat$Sex+Test_Patient_dat$Age)$residuals

plot(Test_Patient_dat$Age,cPAD_rvm,xlim=c(10,75),ylim=c(-40,40),pch=16,
     xlab="Chronological age (Test Patient)",ylab="Brain PAD corrected for age")
abline(lm(cPAD_rvm~Test_Patient_dat$Age))
abline(h=0,lty=2,col="red")
mtext("D", side=3, line=-27, cex=2, adj = 0.5, outer=TRUE)

dev.off()
## png 
##   2

Brain PAD correction -table

##### Brain PAD correction 
test_dat <- rbind(Test_HC_data,Test_Patient_dat)
test_dat$Dx <- c(rep(0,nrow(Test_HC_data)),rep(1,nrow(Test_Patient_dat)))
test_dat$Dx <- factor(test_dat$Dx)

pred_age <- predict(rvm.fit.77 , as.matrix(test_dat[,c(3,5:81)]))
test_error_rvm <- pred_age-test_dat$Age
PAD_rvm <- test_error_rvm
mean(PAD_rvm) ; sd(PAD_rvm) ; median(PAD_rvm)
## [1] -1.815347
## [1] 12.19153
## [1] -0.6313292
Model_linear <- lm(PAD_rvm~test_dat$Age+test_dat$Sex)
NullModel <- lm(PAD_rvm~1)
anova(NullModel,Model_linear,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ test_dat$Age + test_dat$Sex
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    441 65547                           
## 2    439 28222  2     37325 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
test_dat$Age2 <- test_dat$Age^2
Model_quad <- lm(PAD_rvm~test_dat$Sex+test_dat$Age+test_dat$Age2)
anova(NullModel,Model_quad,test="Chisq") # Reject H0 (Full model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rvm ~ 1
## Model 2: PAD_rvm ~ test_dat$Sex + test_dat$Age + test_dat$Age2
##   Res.Df   RSS Df Sum of Sq  Pr(>Chi)    
## 1    441 65547                           
## 2    438 28153  3     37394 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(Model_linear,Model_quad,test="Chisq") # Not Reject H0 (Reduced model is better)
## Analysis of Variance Table
## 
## Model 1: PAD_rvm ~ test_dat$Age + test_dat$Sex
## Model 2: PAD_rvm ~ test_dat$Sex + test_dat$Age + test_dat$Age2
##   Res.Df   RSS Df Sum of Sq Pr(>Chi)
## 1    439 28222                      
## 2    438 28153  1    68.723   0.3011
cPAD_rvm <- lm(PAD_rvm~test_dat$Sex+test_dat$Age)$residuals
mean(cPAD_rvm) ; sd(cPAD_rvm) ; median(cPAD_rvm)
## [1] -4.868136e-16
## [1] 7.999728
## [1] -0.8784845
library(plotrix)
dat <- data.frame(SubjID=test_dat$SubjID,Brain_PAD = c(PAD_rvm),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))

dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"   
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

library(dplyr)
palette <- c("FE-SSD"="darkgreen","TRS"="red","Control"="grey","Schizophrenia"="blue","Patient"="orange")

dat %>% dplyr::select(Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 2 × 4
##   Status    mean    sd std.error
##   <chr>    <dbl> <dbl>     <dbl>
## 1 Control -3.70   12.1     0.834
## 2 Patient -0.128  12.1     0.792
dat %>% dplyr::select(Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_sch     mean    sd std.error
##   <chr>         <dbl> <dbl>     <dbl>
## 1 Control       -3.70  12.1     0.834
## 2 Schizophrenia -1.51  11.6     0.831
## 3 <NA>           7.20  11.9     1.96
dat %>% dplyr::select(Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_trs   mean    sd std.error
##   <chr>       <dbl> <dbl>     <dbl>
## 1 Control    -3.70   12.1     0.834
## 2 TRS        -2.65   11.4     1.69 
## 3 <NA>        0.493  12.2     0.891
dat %>% dplyr::select(Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_ssd  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -3.70  12.1     0.834
## 2 FE-SSD      5.39  12.3     1.45 
## 3 <NA>       -2.59  11.2     0.881
ggplot(dat,aes(x = Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette) +
  coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_sch),],aes(x = Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_ssd),],aes(x = Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_trs),],aes(x = Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

summary(lm(Brain_PAD~Status+test_dat$Sex+test_dat$Age,data=dat))
## 
## Call:
## lm(formula = Brain_PAD ~ Status + test_dat$Sex + test_dat$Age, 
##     data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.113  -5.046  -0.857   4.577  33.873 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   15.61240    1.55032  10.070  < 2e-16 ***
## StatusPatient  4.04615    0.74835   5.407 1.06e-07 ***
## test_dat$Sex   3.73670    0.76567   4.880 1.49e-06 ***
## test_dat$Age  -0.67469    0.02701 -24.978  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.772 on 438 degrees of freedom
## Multiple R-squared:  0.5964, Adjusted R-squared:  0.5936 
## F-statistic: 215.7 on 3 and 438 DF,  p-value: < 2.2e-16
summary(lm(Brain_PAD~status_sch+test_dat$Sex+test_dat$Age,data=dat))
## 
## Call:
## lm(formula = Brain_PAD ~ status_sch + test_dat$Sex + test_dat$Age, 
##     data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.243  -5.090  -0.875   4.716  33.706 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             14.95170    1.59412   9.379  < 2e-16 ***
## status_schSchizophrenia  4.03490    0.78702   5.127 4.59e-07 ***
## test_dat$Sex             3.96796    0.80572   4.925 1.24e-06 ***
## test_dat$Age            -0.66716    0.02858 -23.342  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.74 on 401 degrees of freedom
##   (결측으로 인하여 37개의 관측치가 삭제되었습니다.)
## Multiple R-squared:  0.5797, Adjusted R-squared:  0.5765 
## F-statistic: 184.3 on 3 and 401 DF,  p-value: < 2.2e-16
dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]
##     SubjID Brain_PAD  Status status_ssd status_trs    status_sch
## 415    T19 12.636224 Patient     FE-SSD        TRS Schizophrenia
## 419    T24  4.598935 Patient     FE-SSD        TRS Schizophrenia
dat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA

table(dat_for_regression$status_trs)
## 
## Control     TRS 
##     209      46
table(dat_for_regression$status_ssd)
## 
## Control  FE-SSD 
##     209      72
test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])

dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]

summary(lm(Brain_PAD~dat_for_regression$status_ssd+test_dat_for_regression$Sex+
             test_dat_for_regression$Age,data=dat_for_regression))
## 
## Call:
## lm(formula = Brain_PAD ~ dat_for_regression$status_ssd + test_dat_for_regression$Sex + 
##     test_dat_for_regression$Age, data = dat_for_regression)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.0574  -4.7342  -0.7179   4.3585  28.0182 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         18.55111    1.65831  11.187  < 2e-16 ***
## dat_for_regression$status_ssdFE-SSD  3.30188    1.00980   3.270 0.001192 ** 
## dat_for_regression$status_ssdTRS     5.53278    1.20911   4.576 6.78e-06 ***
## test_dat_for_regression$Sex          3.13507    0.84218   3.723 0.000233 ***
## test_dat_for_regression$Age         -0.72666    0.03006 -24.177  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.18 on 322 degrees of freedom
##   (결측으로 인하여 117개의 관측치가 삭제되었습니다.)
## Multiple R-squared:  0.6766, Adjusted R-squared:  0.6726 
## F-statistic: 168.4 on 4 and 322 DF,  p-value: < 2.2e-16
# ANOVA results (Control vs FE-SSD vs TRS)
summary(res_aov <- aov(dat_for_regression$Brain_PAD ~ dat_for_regression$status_ssd))
##                                Df Sum Sq Mean Sq F value   Pr(>F)    
## dat_for_regression$status_ssd   2   4483  2241.7    15.5 3.72e-07 ***
## Residuals                     324  46850   144.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.
pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "fdr")$p.value
##             Control       FE-SSD
## FE-SSD 2.004933e-07           NA
## TRS    5.943715e-01 0.0006832325
pairwise.t.test(dat_for_regression$Brain_PAD,dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value
##             Control      FE-SSD
## FE-SSD 2.004933e-07          NA
## TRS    1.000000e+00 0.001366465
# t.test results (Control vs Schizo)
t.test(dat$Brain_PAD~dat$status_sch)
## 
##  Welch Two Sample t-test
## 
## data:  dat$Brain_PAD by dat$status_sch
## t = -1.8556, df = 402.68, p-value = 0.06425
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -4.5005093  0.1299218
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -3.696547                   -1.511253
t.test(dat$Brain_PAD~dat$status_sch, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  dat$Brain_PAD by dat$status_sch
## t = -1.8534, df = 403, p-value = 0.06455
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -4.5031844  0.1325969
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -3.696547                   -1.511253
#### corrected_Brain_PAD ####

dat <- data.frame(SubjID=test_dat$SubjID,corrected_Brain_PAD = c(cPAD_rvm),Status=rep(c("Control","Patient"),c(nrow(Test_HC_data),nrow(Test_Patient_dat))))

dat$status_sch <- dat$status_trs <- dat$status_ssd <- rep(NA,nrow(dat))
dat$status_sch[test_dat$SubjID %in% Test_HC_data$SubjID ] <- "Control"   
dat$status_sch[test_dat$SubjID %in% Test_schizophrenia_dat$SubjID ] <- "Schizophrenia"
dat$status_trs[test_dat$SubjID %in% Test_TRS_dat$SubjID ] <- "TRS"
dat$status_trs[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

dat$status_ssd[test_dat$SubjID %in% Test_FESSD_dat$SubjID ] <- "FE-SSD"
dat$status_ssd[test_dat$SubjID %in% Test_HC_data$SubjID  ] <- "Control"

dat %>% dplyr::select(corrected_Brain_PAD,Status) %>% group_by(Status) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 2 × 4
##   Status   mean    sd std.error
##   <chr>   <dbl> <dbl>     <dbl>
## 1 Control -2.09  6.47     0.447
## 2 Patient  1.87  8.76     0.574
dat %>% dplyr::select(corrected_Brain_PAD,status_sch) %>% group_by(status_sch) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_sch     mean    sd std.error
##   <chr>         <dbl> <dbl>     <dbl>
## 1 Control       -2.09  6.47     0.447
## 2 Schizophrenia  1.80  8.89     0.635
## 3 <NA>           2.25  8.15     1.34
dat %>% dplyr::select(corrected_Brain_PAD,status_trs) %>% group_by(status_trs) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_trs  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -2.09  6.47     0.447
## 2 TRS         3.17  9.04     1.33 
## 3 <NA>        1.55  8.68     0.635
dat %>% dplyr::select(corrected_Brain_PAD,status_ssd) %>% group_by(status_ssd) %>% summarise_each(funs(mean,sd,std.error))
## # A tibble: 3 × 4
##   status_ssd  mean    sd std.error
##   <chr>      <dbl> <dbl>     <dbl>
## 1 Control    -2.09  6.47     0.447
## 2 FE-SSD      1.66  7.88     0.929
## 3 <NA>        1.97  9.14     0.721
ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette) +
  coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

library(ggplot2)
ggplot(dat,aes(x = corrected_Brain_PAD, fill = Status)) + geom_density(alpha = 0.5) + scale_fill_manual(values=palette[c(3,5)]) +
  coord_cartesian(xlim =c(-30,30)) + theme_bw()

A <- ggplot(dat[!is.na(dat$status_sch),],aes(x = corrected_Brain_PAD, fill = status_sch)) + scale_fill_manual(values=palette[c(3,4)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

B <- ggplot(dat[!is.na(dat$status_ssd),],aes(x = corrected_Brain_PAD, fill = status_ssd)) + scale_fill_manual(values=palette[c(3,1)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

C <- ggplot(dat[!is.na(dat$status_trs),],aes(x = corrected_Brain_PAD, fill = status_trs)) +scale_fill_manual(values=palette[c(3,2)]) + geom_density(alpha = 0.5) + coord_cartesian(xlim =c(-30,30)) + theme_bw() + guides(fill=guide_legend(title=NULL))

library(ggpubr)
png("S6_2(RVR).png", width = 1400, height = 300)
ggarrange(A,B,C,
          labels = c("A","B","C"), 
          ncol = 3, nrow = 1,
          font.label = list(size=20))
dev.off()
## png 
##   2
# ANOVA results (Control vs FE-SSD vs TRS)
dat_for_regression <- dat
replicated <- setdiff(which(!is.na(dat$status_trs) & !is.na(dat$status_ssd)),1:214)
dat[replicated,]
##     SubjID corrected_Brain_PAD  Status status_ssd status_trs    status_sch
## 415    T19            3.254155 Patient     FE-SSD        TRS Schizophrenia
## 419    T24           -1.427691 Patient     FE-SSD        TRS Schizophrenia
dat_for_regression <- rbind(dat_for_regression,dat[replicated,])
dat_for_regression[replicated,]$status_ssd <- NA
dat_for_regression[(nrow(dat_for_regression)-1):nrow(dat_for_regression),]$status_trs <- NA

table(dat_for_regression$status_trs)
## 
## Control     TRS 
##     209      46
table(dat_for_regression$status_ssd)
## 
## Control  FE-SSD 
##     209      72
test_dat_for_regression <- test_dat
test_dat_for_regression <- rbind(test_dat_for_regression,test_dat[replicated,])

dat_for_regression$status_ssd[is.na(dat_for_regression$status_ssd)] = dat_for_regression$status_trs[is.na(dat_for_regression$status_ssd)]

summary(res_aov <- aov(dat_for_regression$corrected_Brain_PAD~
                         dat_for_regression$status_ssd))
##                                Df Sum Sq Mean Sq F value   Pr(>F)    
## dat_for_regression$status_ssd   2   1483   741.3   14.31 1.11e-06 ***
## Residuals                     324  16786    51.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 117개의 관측치가 삭제되었습니다.
pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
               dat_for_regression$status_ssd, p.adj = "fdr")$p.value
##             Control    FE-SSD
## FE-SSD 2.479033e-04        NA
## TRS    3.036364e-05 0.2675513
pairwise.t.test(dat_for_regression$corrected_Brain_PAD,
               dat_for_regression$status_ssd, p.adj = "bonferroni")$p.value
##             Control   FE-SSD
## FE-SSD 4.958066e-04       NA
## TRS    3.036364e-05 0.802654
# t.test results (Control vs Schizo)
t.test(dat$corrected_Brain_PAD~dat$status_sch)
## 
##  Welch Two Sample t-test
## 
## data:  dat$corrected_Brain_PAD by dat$status_sch
## t = -5.0093, df = 354.77, p-value = 8.63e-07
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -5.416722 -2.362566
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -2.088031                    1.801613
t.test(dat$corrected_Brain_PAD~dat$status_sch, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  dat$corrected_Brain_PAD by dat$status_sch
## t = -5.0592, df = 403, p-value = 6.411e-07
## alternative hypothesis: true difference in means between group Control and group Schizophrenia is not equal to 0
## 95 percent confidence interval:
##  -5.401066 -2.378221
## sample estimates:
##       mean in group Control mean in group Schizophrenia 
##                   -2.088031                    1.801613
summary(res_aov <- aov(dat$corrected_Brain_PAD~dat$status_sch))
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## dat$status_sch   1   1530  1530.3   25.59 6.41e-07 ***
## Residuals      403  24094    59.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 결측으로 인하여 37개의 관측치가 삭제되었습니다.